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24th  Scandinavian  Symposium  on  Physical  Acoustics 
Ustatoset  28  January  -  31  January  2001 

Preface 


The  Scandinavian  Symposium  on 
Physical  Acoustics  was  held  at 
Ustaoset,  Norway  in  the  period  28th  of 
January  to  the  31st  of  January  2001. 

The  meeting  represented  the  24th  in  a 
series  of  Scandinavian  co-operation 
meetings  in  acoustics/hydrodynamics, 
and  was  arranged  by  the  acoustics  and 
optics  group  of  the  Norwegian 
Physical  Society.  The  symposium  was 
financially  supported  by  the  Office  of 
Naval  Research  International  Field 
Office  (Navy  Grant  N00014-01-1- 
1007).  For  this  we  are  most  grateful.  * 

52  participants  from  Canada, 

Denmark,  England,  Lithuania,  Sweden, 
the  United  States  of  America,  and 
Norway  attended  the  meeting.  A  total 
of  27  papers  were  presented.  They 
covered  a  wide  range  of  topics: 
underwater  acoustics,  new  solutions  to 
sound  propagation  problems, 
interaction  of  sound  and  porous  media, 
non-linear  acoustics,  transducer 
technology,  and  applications  of 
ultrasound  in  industry.  The  papers 
were  presented  in  5  oral  sessions.  The 
presentation  time  allowed  for  each 
presentation  was  25  minutes. 

A  tradition  in  these  symposia  has  been 
to  leave  some  time  around  noon  to  go 
skiing,  so  also  this  year.  The  weather 
was  good,  and  the  nearby  “lpyper”  and 
downhill  slopes  at  Geilo  offered 
possibilities  at  all  levels  of  ambition. 
This,  followed  by  lunch,  allowed  the 


participants  to  meet  fresh  for  the 
afternoon  sessions. 

It  is  the  aim  of  the  organisers  to  attract 
students  at  master  and  doctoral  levels 
to  present  their  work.  The  Office  of 
Naval  Research  International  Field 
Office  and  The  Norwegian  Physical 
Society  must  be  acknowledged  for 
their  financial  support  to  students  who 
attended  the  meeting. 

The  proceedings  contain  18  of  the 
presentations.  They  appear  in  the  order 
they  were  presented  during  the 
symposium.  A  separate  index  for  the 
papers  has  been  added.  There  was  no 
page  limitation  for  the  manuscripts. 

The  organisers  would  like  to  thank 
sincerely  the  speakers  and  participants 
for  making  the  meeting  an  interesting 
and  enjoyable  event.  Special  thanks  to 
the  session  chairpersons,  who 
prevented  the  sessions  to  stray  too  far 
from  the  time  schedule. 

The  Acoustics  and  Optics  Group  of  the 
Norwegian  Physical  Society  plans  to 
organise  a  new  meeting  in  the  year 
2002.  This  will  be  the  25th  in  the 
series!  Invitations  will  be  sent  early  in 
the  autumn  of  2001. 

Ulf  Kristiansen  was  the  co-ordinator  of 
the  Ustaoset  2001  meeting. 

The  content  of  this  document  does  not 
necessarily  reflect  the  position  or  the  policy  of 
the  United  States  government. 
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department  of  Physics,  University  of  Bergen,  Norway 
2  Nansen  Environmental  and  Remote  Sensing  Center 
Edv.  Griegsv.  3A,  N-5059  Bergen,  Norway 
e-mail:  halvor.hobak@fi.uib.no 


Introduction 

The  objective  of  the  research  project  “Acoustic  Monitoring  of  the  Ocean  Climate  in  the  Arctic”  (AMOC)[l] 
is  to  investigate  the  possibility  of  using  acoustic  methods  to  monitor  temperature  and  ice  thickness  in  the 
Arctic  Ocean  on  a  long  term  basis,  with  great  relevance  to  global  warming.  Since  main  influx/efflux  of 
water  (and  heat)  to  the  Arctic  Ocean  is  through  the  Fram  Strait,  monitoring  the  temperature  distribution 
and  water  flux  in  the  Fram  Strait  is  an  important  issue  in  this  respect.  In  this  presentation  we  focus 
on  certain  aspects  of  acoustic  thermometry  using  ray  tracing.  In  order  to  apply  acoustic  means  for 
monitoring  temperature  and  flux  of  water  masses  it  is  essential  to  be  able  to  identify  arrival  of  individual 
rays  at  the  receiver  location,  while  the  environment  undergoes  seasonal  and  shorter  time  scale  variations. 
The  present  study  is  based  on  ray  simulations  using  environmental  data  for  the  period  1950-1990  in  terms 
of  decadal  means  separated  in  winter  (January  -  March)  and  summer  (May-July)  seasons.  In  addition  a 
detailed  data  set  of  one  passing  meso  scale  eddy  was  used.  The  method  applied  is  analysis  of  "ray  fans" 
obtained  by  launching  a  large  number  of  rays  in  a  narrow  fan  centered  about  the  horizontal.  The  main 
conclusion  is  that  some  rays  may  be  identified  most  of  the  time,  even  during  passing  eddies. 


Fram  Strait  Environment 

The  environmental  data  used  in  this  project  is  prepared  from  the  Joint  US-Russian  Arctic  Oceanography 
Atlas  CDs  [2]  and  is  gridded  for  the  Fram  Strait  in  9  stations  separated  by  about  2.8  degrees  along  the 
79°  N  latitude,  covering  the  longitude  11°  W  to  11°  E.  About  200  km  of  the  600  km  wide  Fram  Strait 
between  Greenland  and  Spitzbergen  at  79°  N  is  more  than  2000m  deep.  The  bathymetry  of  the  strait 
is  shown  in  Fig.  1,  with  an  interpolated  color  map  of  the  mean  sound  speed  for  the  summer  seasons 
1950-1989.  The  width  of  the  strait  allows  influx  of  warm  water  to  the  Arctic  Ocean  and  outflux  of  cold 
water  and  ice  to  be  separated  horizontally,  with  the  warm  water  flowing  northwards  on  the  Spitzbergen 
side  and  cold  water  flowing  southwards  on  the  Greenland  side,  which  is  ice-  covered  almost  all  year  round. 
This  results  in  strong  horizontal  sound  speed  gradients  in  addition  to  vertical,  as  clearly  seen  in  Fig.  1. 
Here  Greenland  is  at  the  left  and  Spitzbergen  at  the  right  hand  side,  and  the  origin  is  placed  at  11°  W. 
The  seasonal  variations  occur  mostly  in  the  upper  600  m  of  the  water  column.  From  an  acoustical  point 
of  view  the  Fram  Strait  is  an  extremely  difficult  environment  to  monitor  because  of  the  two  opposing 
currents  which  is  displaced  horizontally.  In  this  study  we  have  therefore  focussed  on  placing  the  sound 
source  in  the  middle  of  the  strait,  and  a  receiver  array  on  the  the  eastern  flank  near  Spitzbergen.  In  the 
examples  shown  the  source  is  located  at  250  km,  and  the  receiver  at  410  km.  Because  of  the  the  sound 
ducts  which  appear  in  the  upper  layers  it  is  necessary  to  lower  the  source  to  300  -  500m  depth  -  in  the 
examples  shown  500m.  The  averaged  environmental  data  contain  no  information  of  short  term  variations 
like  passing  meso  scale  eddies  (time  scale  weeks),  tidal  currents  or  internal  waves.  In  order  to  model  a 
passing  eddy  data  from  the  AOGC  expedition  in  1997  [3]  has  been  used  in  combination  with  the  averaged 
environmental  data. 
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Figure  1:  Bathymetry  and  summer  mean  sound  speed  map  of  the  Pram  Strait 


Ray  fans 

A  preliminary  study  was  made  by  using  the  ray  tracing  code  "RAY"  [4]  to  search  for  "eigenrays"  between 
the  source  and  a  receiver  at  300m  depth.  For  the  different  seasons  and  decades  the  number  of  eigenrays 
varied  between  2  and  14.  A  study  of  the  eigenray  signatures,  determined  by  launch  angle,  number  of 
top  and  bottom  bounces  and  angle  at  receiver,  revealed  a  marked  difference  between  summer  and  winter 
seasons,  and  a  rather  confusing  pattern  emerged.  In  order  to  investigate  this  further  a  different  approach 
was  taken.  Instead  of  searching  for  eigenrays  a  large  number  of  rays  (1200)  was  launched  within  a  narrow 
angle  interval  about  the  horizontal  (ray  fan),  and  the  depth  and  travel  time  at  the  receiver  range  was 
recorded.  Only  rays  actually  reaching  the  receiver  range  was  kept.  A  plot  of  the  ray  fan  for  summer  50s 
is  shown  in  Fig.  3  as  ray  depth  at  receiver  against  source  angle  (negative  angles  point  downwards),  Fig. 
4  as  travel  time  against  source  angle  and  Fig.  5  as  depth  against  arrival  time  (these  are  3  projections  of 
a  3-D  curve  in  the  space:  source  angle,  ray  depth  and  travel  time  -  a  3-D  plot  is  interesting,  but  does  not 
really  reveal  important  new  information). 

Each  ray  arrival  is  plotted  as  a  circle,  and  the  density  of  these  represent  the  intensity  of  the  sound. 
The  horizontal  line  indicates  a  receiver  at  300m  depth.  It’s  intersections  with  the  curve  identifies  the 
eigenrays  for  this  receiver.  The  eigenray  algorithm  in  RAY  found  5  eigenrays,  corresponding  to  the 
densest  intersections.  Do  observe  that  in  the  "depth  against  time"  plot,  Fig.  3,  the  deepest  rays  are 
the  last  to  arrive!  This  is  found  to  be  the  case  in  general  (with  only  one  exception),  and  is  contrary 
to  ordinary  sound  channel  behavior.  The  fastest  rays  -  in  most  cases  a  narrow  ray  fan  launched  almost 
horizontally  (in  Fig.  3  -1°  to  +0.5  °)  -  propagate  in  a  narrow,  undulating  band  about  an  almost  straight 
line  from  the  source  to  the  receiver,  which  is  reached  at  depths  varying  from  800  m  to  50  m.  The  case 
summer  50  is  shown  in  Fig.  5.  Thus,  these  rays  seem  to  follow  the  range  dependent  sound  channel  axis. 
The  reason  why  they  are  the  first  to  arrive  must  be  due  to  the  path  length  being  as  much  shorter  than 
for  the  deeper  rays  that  it  compensate  for  the  lower  sound  speed. 

Ray  fans  plotted  as  in  Fig.  2-4  are  helpful  for  investigating  stability  of  eigenrays.  As  seasons  and 
decades  vary  the  patterns  are  deformed  and  shifted  -  a  really  bad  case  of  such  deformation  is  shown  in 
Fig.  6,  for  winter  60s.  Here  the  eigenray  algorithm  in  RAY  found  only  2  eigenrays  (it  turns  out  that  in 
this  case  the  ray  fan  looks  more  normal  if  the  source  is  shifted  to  shallower  depths,  300-400m,  but  then 
the  deepest  rays  penetrate  only  down  to  600m). 

In  an  attempt  to  obtain  a  stable  reference  frame  for  the  ray  fan  results  Fig.  7  shows  a  "generic  plot" 
of  this  type,  obtained  by  disregarding  range  dependence  (using  the  sound  speed  profile  at  the  source)  for 
summer  60. 

Evidently,  the  ray  fan  is  rather  symmetric  about  zero  source  angle,  and  shows  3  undulations  to  each 
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Figure  5:  Trace  of  fastest  rays,  summer  50 


side.  These  are  narrowly  peaked  at  the  top  and  wider  at  the  deep  end.  The  outer  flanks  fall  off  steeply 
towards  depths  about  1  km.  As  seasons  vary  both  amplitudes  and  widths  of  these  undulations  fluctuate, 
in  particular  in  the  region  ±2°,  and  sometimes  they  are  severely  deformed  (in  Fig.  2  the  rightmost 
undulation  is  actually  No  2,  not  3).  The  pattern  may  be  shifted  upwards  or  downwards,  and  some 
undulations  may  merge  (in  particular  the  central  ones)  or  split.  The  trend  emerging  from  such  analysis 
is  that  the  rays  belonging  to  the  outer  flanks  are  the  ones  most  likely  to  be  recognizable  at  all  times, 
while  those  launched  more  horizontally  are  less  stable.  Also,  since  the  stablest  rays  are  also  the  deepest, 
they  are  the  most  promising  candidates  for  acoustical  thermometry  in  the  FVam  Strait.  For  convenience 
these  eigenrays  are  labeled  "a,b"and  "g,h",  where  a  and  h  are  on  the  outermost  flanks,  and  b  and  g  on 
the  inner  side  of  the  outer  undulations  (in  Fig.  2  only  eigenray  "b"  is  present,  near  -4°  source  angle).  In 
Fig.  8  is  shown  traces  of  these  eigenrays  for  all  seasons  (source  depth  500m,  receiver  depth  300m).  They 
cover  the  water  column  down  to  1  km  and  more,  and  thus  the  dominating  part  of  the  West  Spits-bergen 
Current. 

The  ray  fan  plots  also  show  that  one  should  have  receivers  at  several  depths,  in  order  to  increase  the 
probability  of  finding  stable  eigenrays  -  if  they  fail  at  one  depth  for  a  period  there  would  be  some  at  other 
depths,  and  vice  versa. 

In  order  to  demonstrate  the  sensitivity  of  acoustic  thermometry  based  on  the  stable  rays  the  travel 
times  of  these  rays  are  shown  in  Fig.  9.  There  is  a  marked  difference  between  summer  and  winter  seasons; 
typically  about  100  ms  faster  in  the  summer.  Therefore  the  seasons  are  separated  in  the  figure.  Note  that 
ray  h  is  missing  in  all  summer  seasons.  The  eigenrays  show  basically  the  same  trend:  Slightly  decreasing 
travel  times  with  time,  and  more  pronounced  at  the  later  decades.  The  change  is  small,  however:  Winter 
rays  arrive  31±2ms  earlier  in  1980s  than  in  1959s,  corresponding  to  an  increase  in  sound  speed  0.42  ±0.03 
m/s,  or  an  increase  in  temperature  of  roughly  0.09°C.  Inversion  based  on  details  of  the  ray  paths  has 
not  yet  been  performed. 


Meso  scale  eddies 

The  original  sound  speed  data  for  one  meso  scale  eddy  taken  by  the  1997  AOGC  expedition  [3]  is  shown 
in  Fig.  10.  This  data  set  was  incorporated  into  the  environmental  data  in  the  following  way.  First, 
sound  speed  data  was  interpolated  on  a  regular  grid  horizontally  and  vertically  in  terms  of  sound  speed 
differences  (global  mean  =  zero).  Secondly,  margins  were  added  to  both  sides  and  bottom  with  tapered 
sound  speed  differences,  in  order  to  avoid  artifacts  due  to  discontinuities.  Finally  the  data  set  was  merged 
onto  the  environmental  data  sets  (added),  with  an  amplitude  factor  adjustable  between  0  and  1.  In  this 
way  one  could  to  some  extent  simulate  a  passing  eddy. 

As  the  meso  scale  eddy  is  increased  the  ray  fans  become  distorted,  first  near  zero  source  angle,  later 
also  on  the  outer  undulations.  An  example  is  shown  in  Fig.  11  for  maximum  meso  amplitude,  summer 
60.  The  pattern  tends  to  become  chaotic,  and  the  number  of  eigenrays  can  become  very  high,  but  with 
forbiddingly  low  amplitude.  However,  the  "stable"  eigenrays,  a,b,g  and  h  turn  out  to  be  less  influenced 
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Winter  60s,  source  250  km,  500m,  range  410  km 


Figure  10:  Map  of  the  ssp-data  for  the  meso  scale  eddy. 


than  the  rest,  and  are  therefore  promising  candidates  for  eigenrays  for  acoustical  thermometry  in  a  future 
experiment. 


S60,mcso  factor  1.0,  Source  250  km,  range  410  km.  Depth  vs  arrival  time 


Source  angle  -  deg 

Figure  11:  Ray  fan  meso  factor  1,  summer  60,  depth  against  source  angle 


Conclusions 

Ray  fan  analysis  is  a  valuable  tool  for  assessing  eigenray  stability.  The  steepest,  deepest  going  rays  are 
the  most  promising  candidates  for  acoustic  thermometry  in  the  Pram  Strait  over  a  prolonged  period. 
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Abstract 

A  field  trial  has  been  conducted  in  the  Stockholm  Archipelago  in  the  Baltic  Sea  with  the  objective  of  investigating  the  influence  of 
the  waveform  on  the  reverberation  level.  The  performance  of  three  LFM  pulses  of  different  frequency  content  and  bandwidth  and 
one  20  kHz  CW  pulse  was  compared  after  matched  filter  processing.  The  more  broadband  signals  have,  as  expected  from  theory,  a 
lower  reverberation  level.  A  simulation  study  was  also  performed  and  the  agreement  between  experimental  data  and  the  model  is 
satisfactory  for  ranges  up  to  a  few  hundred  meters,  above  which  the  received  levels  were  below  the  noise  floor  of  the  experiment 


Shallow  water  reverberation 

Active  sonar  performance  is  often  reverberation  limited  in 
shallow  waters.  In  the  Baltic  Sea  a  common  bottom  type 
consists  of  bedrock  covered  by  soft  sediments  of  varying 
thickness.  The  sediment  can  be  penetrated  by  signals, 
especially  of  low  frequency,  and  be  reflected  either  in 
sediment  layer  boundaries  or  in  the  sediment-bedrock 
interface.  The  limited  depth  of  the  Baltic  Sea  enhances  the 
reverberation  level  experienced  in  active  sonar  applications. 
To  maximize  the  performance  of  active  sonar  systems, 
relatively  high  frequencies  have  traditionally  been  chosen 
for  Baltic  Sea  applications. 

In  general,  there  are  several  parameters  that  have  impact 
on  the  reverberation  level.  Some  are  due  to  the 
environment.  This  includes  the  bottom  topography  and  the 
bottom  type  of  the  site,  and  the  sound  speed  profile.  There 
are  also  sonar  parameters  that  are  important,  such  as  the 
lobe  width.  This  work  is  a  study  of  how  the  waveform, 
especially  the  pulse  length,  influences  the  reverberation 
level.  Since  we  use  a  matched  filter  when  processing  our 
data,  the  pulse  length  is  determined  by  the  signal  bandwidth. 
A  field  trial  was  conducted  in  a  shallow  bay  in  the 
Stockholm  Archipelago  and  the  measured  reverberation 
levels  from  a  few  pulse  types  have  been  compared  with  a 
theoretical  model.  More  details  about  this  work  is  found 

in  rn- 

Field  trial 

The  field  trial  was  conducted  in  May  2000.  The 
experiments  were  carried  out  at  FOI’s  permanent  field 
laboratory  situated  in  the  southern  part  of  the  Stockholm 
Archipelago.  At  this  site  there  is  a  floating  platform  from 
which  the  transmitter  was  deployed.  Around  the  transmitter 
a  ring  with  six  receiver  hydrophones  was  mounted.  This 
means  that  the  measurements  were  carried  out  in  monostatic 


mode.  The  transmitter,  a  TOPAS  120  from  Simrad,  has  a 
narrow  lobe  width,  typically  4  degrees  at  the  frequencies 
used.  The  lobe  width  has  a  weak  frequency  dependency, 
and  is  narrower  for  higher  frequencies.  The  source  level  of 
this  parametric  sonar  is  rather  low,  around  1 90  dB  re  1 
pPa  at  1  m  distance  at  the  frequencies  used.  In  this 
frequency  range  the  ring  of  receiver  hydrophones  has  an 
almost  omm-directional  response.  For  more  details  on  the 
equipment,  see  [l]-[2] . 

The  transmitter  and  receiver  hydrophones  were  kept  at  a 
depth  of  1 .3  m  throughout  the  measurements.  The  tilt  angle 
of  the  transmitter  and  receiver  was  varied  from  1  degree 
upwards  to  10  degrees  downwards.  Horizontally,  a  sector 
spanning  60  degrees  was  covered.  This  sector  contains 
some  variation  in  the  bottom  topography,  and  the  influence 
of  this  was  probed.  The  bottom  type  can  be  considered 
constant  in  the  field  trial  area,  and  it  consists  of  a  smooth 
clay-covered  seabed. 

In  active  sonar  applications  it  is  the  received  signal  to 
reverberation  level  that  decides  the  performance.  To  have 
known  target  echoes  to  compare  the  reverberation  levels 
with,  two  comer  reflectors  of  0.5  m  dimension  were 
deployed.  One  was  at  a  distance  of  100  m  from  the 
transmitter/receiver,  while  the  second  one  was  placed  at 
about  330  m  range.  The  echo  from  the  latter  reflector  was 
in  most  cases  too  weak  to  be  observed. 

The  measurements  were  carried  out  during  three  days  with 
some  variations  in  the  sound  speed  profile.  While  the  profile 
was  downward  refracting  during  the  first  two  measurement 
days,  it  changed  into  iso-velocity  the  last.  The  day  to  day 
variations  in  the  sound  speed  profile  had  a  larger  impact 
on  the  measured  reverberation  levels  than  the  ones  due  to 
variations  in  the  horizontal  angle,  and  thus  the  seabed 
topography. 
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Figure  1:  Received  reverberation  levels  for  the  four  pulse-types 
used  in  this  study.  The  two-way  travel  time  has  been  translated 
into  an  equivalent  distance  from  the  transmitter/receiver.  This 
figure  shows  the  combined  result  from  all  tilt  angles  when  the 
transmitter/receiver  was  aimed  towards  the  closer  comer  reflector, 
in  total  about  200  pings  per  pulse. 


Figure  2:  A  close  view  of  the  echo  from  a  comer  reflector 
deployed  to  study  the  echo  to  reverberation  level  for  the  different 
waveforms. 


Analysis  of  the  data 

The  expected  bottom  reverberation  level  is  proportional 
to  the  bottom  area  insonified  by  the  emitted  pulse  [3].  This 
area  depends  on  the  beam  width,  the  range,  the  sound  speed 
and  the  pulse  length.  There  is  thus  a  proportionality  between 
the  received  reverberation  level  and  the  pulse  length,  r. 
The  recorded  data  set  was  analysed  using  a  matched  filter, 
i.e.  a  synthetic  replica  of  the  transmitted  signal  was 
correlated  with  the  received  echoes.  This  process  means  a 
pulse  compression  in  time  with  a  factor  B  r,  where  B  is  the 
bandwidth  of  the  emitted  pulse.  The  result  is  a  decrease  of 
the  reverberation  level  with  a  factor  10-log10(B  r).  In  other 
words,  the  expectation  is  that  a  wider  bandwidth  will  have 
a  shorter  effective  pulse  length  after  the  matched  filter,  and 
this  will  decrease  the  reverberation  level.  To  investigate 
this  hypothesis  empirically,  we  used  a  few  different  pulses: 
one  monofrequency  (CW)  20  kHz  signal,  and  three  linear 
frequency  modulated  (LFM)  pulses  of  varying  frequencies, 
4-8  kHz,  16-20  kHz  and  19-20  kHz.  The  result  of  the 
measurements  is  summarized  in  Figure  1 . 

From  this  figure,  the  expected  behaviour  is  verified.  The 
broadband  signals  have  less  reverberation  than  the 
narrowband-  and  monofrequency  pulses.  The  difference 
between  the  pulse  with  a  bandwidth  of  1  kHz  and  the  ones 
with  4  kHz  bandwidth  is  about  the  expected  value,  6  dB. 

The  early  part  of  the  received  pulse  is  not  shown,  since  the 
receiving  hydrophones  were  mounted  immediately  around 
the  parametric  sonar,  and  were  saturated  while  the  pulse 
was  transmitted.  The  omitted  part  is  slightly  longer  than 
the  pulse  length  to  avoid  erroneous  results. 

The  echo  from  the  close  comer  reflector  is  clearly  visible 
inFigure  1.  This  echo  is  also  shown inFigure  2. 


The  best  echo  to  reverberation  ratio  is  achieved  with  the 
16-20  kHz  LFM  pulse,  with  its  narrow  pulse  width  after 
the  matched  filter.  The  4-8  kHz  LFM  pulse  does  not  resolve 
the  target.  This  is  because  the  wavelength  of  that  pulse  is 
of  the  same  order  as  the  dimension  of  the  target. 

Comparison  with  a  theoretical 
model 

A  two-dimensional  ray-based  model,  SLOPERVB,  has 
been  used  for  this  study.  It  is  based  on  the  model  B1RVB 
[4],  with  the  extension  that  range-dependent  bottom 
topographies  can  be  studied.  The  model  is  described  in 
[1].  The  depth  profiles  used  in  the  simulations  varied  with 
the  bearing  angle  of  the  transmitter/receiver.  Other 
important  environmental  parameters  that  were  adjusted  to 
the  actual  experimental  situation  are  the  sound  speed  profile 
(range  independent)  and  the  bottom  impedance  (a  value  of 
1.9- 10s  kg/m3  was  found  to  be  suitable). 

The  backscattering  at  the  bottom  was  modelled  following 
Lambert’s  law, 

SS  =  10  ■  log10  /I  +10  ■  log10  (sin  Q1  ■  sin  02 ) 

where  SS  is  the  scattering  strength  when  the  grazing 
angles  for  incidence  and  scattering  are  0j  and  02  .The 
value  for  to  the  scattering  strength  that  was  found  to  produce 
the  best  agreement  with  the  measured  data  was 

10  ■  log.0  /J  —  -17dB 

It  is  also  in  agreement  with  other  more  detailed  studies  of 
the  backscattering  strength  at  the  site  [5], 

The  sonar  has  a  lobe  width  with  some  frequency 
dependency,  and  this  was  considered  in  the  computer 
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according  to  the  model  for  three  tilt  angles  of  the  transmitter/ 
receiver  ring. 

modelling.  Since  the  results  from  the  model  are  compared 
with  data  after  a  matched  filter,  the  pulse  length  was 
determined  by  the  bandwidth  of  the  modelled  pulses. 

Figure  3  shows  some  modelling  results.  The  pulse  is  a  1 6- 
20  kHz  LFM  pulse  and  the  tilt  angle  of  the  transmitter/ 
receiver  has  been  varied.  The  result  can  be  compared  with 
Figure  4  which  shows  experimental  data  under  similar 
circumstances.  The  agreement  is  satisfactory,  and  there  are 
some  common  features  in  the  curves.  Please  note  that  the 
experimental  data  contain  an  echo  from  a  comer  reflector, 
and  this  was  not  modelled.  For  larger  ranges,  Le .  late  arrival 
times,  the  model  predicts  lower  levels  than  the  noise  floor 
in  the  experiment,  and  no  comparison  between  experiment 
and  model  can  be  made.  The  absolute  levels  are  also 
compatible  with  each  other. 


Distance  [m] 

Figure  4 :  Similar  to  Figure  3  but  with  experimental  data.  Please 
note  that  the  scales  are  different  in  Figure  3  and  Figure  4. 
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Conclusions 

A  field  trial  designed  to  study  how  reverberation  depends 
on  signal  bandwidth  has  been  conducted.  The  data  have 
been  analysed  using  a  matched  filter  and  the  results  are 
consistent  with  the  theoretical  expectations;  LFM  pulses 
with  4  kHz  bandwidth  give  less  reverberation  than  the  1 
kHz  bandwidth  case.  Data  has  also  been  compared  with  a 
model,  SLOPERVB,  and  for  ranges  up  to  a  few  hundred 
meters  the  agreement  is  satisfactory.  For  longer  ranges  no 
comparison  was  possible  to  do,  since  with  the  present 
source  the  received  reverberation  levels  were  below  the 
noise  floor  of  the  experiment.  The  agreement  between  the 
model  and  the  data  is  encouraging  and  opens  up  the 
possibility  to  do  simulation  studies  and  comparisons  of 
reverberation  levels  from  other  types  of  broadband  signals 
in  the  future,  such  as  binary  phase  shift  keying  (BPSK) 
signals  as  an  example. 
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Abstract 

The  results  of  a  numerical  investigation  of  the  acoustic  propagation  in  the  ice-covered 
region  of  Arctic  Sea  where  the  environment  can  be  considered  range-independent  are 
presented.  Environment  parameters  are  motivated  from  measurements  from  Greenland 
Sea  during  Marginal  Ice  Zone  Experiment,  MIZEX  [4]. 


1  Introduction 

Numerical  modeling  of  acoustic  propagation  in  the  ocean  has  been  in  continuous  development 
throughout  the  last  two  decades.  Before  the  early  1970,  ray  theory  was  applied  for  solving 
propagation  problems  in  underwater  acoustics,  which  is  a  computationally  efficient  but  based 
on  high  frequency  asymptotic  approximation.  Since  early  1970  there  has  been  an  enormous 
expansion  in  computer  science  and  this  progress  has  stimulated  acousticians  to  develop  fre¬ 
quency  dependent  solutions  for  wave  equation.  In  these  numerically  based  techniques  wave 
equation  is  solved  by  normal  mode,  fast  field,  and  parabolic  equation  methods  [3]  .  This  has 
made  it  possible  to  compute  the  acoustic  field  in  complex  ocean  environments. 

Hydrological  measurements  in  the  Marginal  Ice  Zone  (MIZ)  of  the  Greenland  sea  reveals 
the  need  of  range-dependent  models  in  the  ice  edge  zone,  whereas  a  range-independent  hor¬ 
izontal  layered  model  can  be  used  in  ice  covered  regions.  In  this  paper,  the  focus  is  on  the 
ice-covered  region  of  the  Arctic.  SAFARI  numerical  code  [2]  provides  accurate  transmission 
loss  values  for  propagation  in  range-independent  multi-layered  horizontally  stratified  environ¬ 
ments  and  was  therefore  a  useful  tool  for  this  investigation. 

First  solution  of  wave  equation  for  a  range-independent  environment  is  presented  and  mathe¬ 
matical  equations  for  Green’s  function  and  normal  stress  are  given  for  a  3  layered  environment 
(vacuum-solid-fluid-fluid  halfspace)  in  Sec.  2  and  a  numerical  example  (case  1)  is  presented 
to  study  a  deep  ice-covered  sea,  3  other  numerical  examples  study  more  complicated  environ¬ 
ments  in  Sec.  3.  SAFARI  was  applied  for  simulations.  A  summary  ends  the  discussion  in 
Sec.  4  in  which  some  comments  are  given  about  SAFARI. 
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2  Mathematical  basic  equations  for  a  3  layered  model 

By  employing  forward  Fourier  and  Hankel  transform  respectively  on  the  linear  wave  equation 
in  a  range  independent  environment  ,  the  depth-separated,  wave  equation  can  be  obtained  [2]: 

l^-(k,-kl(z)))G(,k,z)  =  ?Q.  (1) 

Where  km  is  medium  wavenumber  and  F(z)  is  forcing  term.  The  total  solution  of  the  depth- 
separated  equation  can  be  expressed  as  a  superposition  of  a  homogeneous  solution  Gjj,  satis¬ 
fying  boundary  conditions,  and  a  particular  solution  Gp,  the  free  field  produced  by  a  source 
in  an  unbounded  medium.  Gp(k,z )  term  is  added  to  the  homogeneous  solution  in  layers 
containing  one  or  more  sources: 


G(k,  z)  =  Gp{k,  z)  +  Gp(k,  z).  (2) 

GH{k,z)  is  found  by  using  the  boundary  conditions.  To  find  the  total  field  at  any  range  r 
and  the  time  response,  we  apply  inverse  Hankel  transform. 


9(r,z)  =  / 
Jo 


G(k,z)Jo(kr)kdk, 


where  J0  is  the  Bessel  function  and  k  is  the  horizontal  wavenumber. 
The  normal  stress  tzz  is  then  given  by: 


tzz  =  -pv2g{r,  z) 


(4) 


where  p  the  fluid  density  ,  u>  the  radial  frequency  and  g(r,  z)  is  the  potential  field. 


In  this  paper  G(k,z)  and  tzz  are  studied.  To  investigate  the  acoustic  propagation  in  the 
ice-covered  region  a  3  layered  model  is  considered  (vacuum,  solid,  fluid,  fluid  half-space). 


The  integral  representation  of  the  homogeneous  solutions  in  the  layer  1  can  be  expressed 
as  :  roo 

9n(r,z)=  [ Ai(k)e-QlZ  +  At(k)ea'z}J0{kr)kdk ,  (5) 

jo 

where  A j  (k)  and  Af  ( k )  are  the  amplitudes  of  the  upgoing  and  downgoing  compressional  field 
in  the  layer  1  and 


<*i 


-y/k2-kl 
iJk2  —  k2, 


k  >  fci, 
k  <  k\. 


roo 

9si(r,z)=  /  [B{(k)e~0lZ  +  B+(k)e^z]J0(kr)kdk, 

JO 


where  = 


( 


i\fk 2  -  k'i, 


k  >  kS) 
k  <  ks, 


B\  (ft)  and  B+ (ft)  are  the  amplitudes  of  the  upgoing  and  downgoing  shear  field  in  the  layer  1, 
9li{riz)  and  gsi(r,z)  are  the  potential  fields  in  the  solid  corresponding  to  the  compressional 
velocity  and  the  shear  velocity  respectively. 
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Using  the  Sommerfeld  radiation  condition,  the  total  solution  in  the  layer  2,  where  the  source 
is  located,  can  be  expressed  in  the  form  of: 

«oo 

g2(r,z)  =  I  {[Az{k)e-a^z-V  +  A+ {k)ea^z-D^}  +  Gp2{k,z,zs)}J0(kr)kdk,  (7) 

Jo 


where  a2  = 


{ 


-yjk2  -  kl  k  >  k2, 
iyjk2  —  k2,  k  <  k2. 


A2  {k)  and  At{k)  are  the  amplitudes  of  the  upgoing  and  downgoing  fields  in  the  layer  2, 
Gp2{k,z,zs)  =  Su  is  the  source  strength  and  D  the  thickness  of  the  layer  1. 

Similarly  the  integral  representation  of  the  solution  in  layer  3  can  be  introduced  as: 


g3{r,z)  =  A+ (k)ea^D+L^J0(kr)kdk, 

Jo 


«3  = 


{ 


—yjk2  —  kl,  k  >  k3, 
iyjk2  -  k^,  k  <  k3, 


(8) 


where  At  is  the  amplitude  of  the  downgoing  field  in  layer  3  and  L  thickness  of  layer  2.  The 
boundary  conditions  for  such  a  model  are  presented  in  [2],  p.  17. 

These  lead  to  a  linear  equation  system  which  can  be  written  as  matrix  form: 


C 


^r 

M 

BT 

B{ 

a2 

At 

At 


—  M, 


where  matrices  C  and  M  are: 


2fcQ!l 

—2ka\ 

-K 

-K 

0 

0 

0 

K 

K 

—2kfii 

- 

2kf31 

0 

0 

0 

piKe~aiD 

P!KeaiD 

—2p\kfJ\e~&xD 

—2pikf3\e^lD 

P2 

P2 

0 

2kaie~aiD 

-2kaieQlD 

-Ke-P'D 

-Ke0lD 

0 

0 

0 

—a\e~aiD 

ke~^D 

ke0xD 

012 

-a2 

0 

0 

0 

0 

0 

p2eot2L 

p2e"2L 

~P3 

0 

0 

0 

0 

—a2e  0,2  L  a2ea2L 

-03 

r 

0 

0 

p2eaAD-**) 

where  K 

=  2k2  -  k2, 

M=- 

S<jj 

4na2 

0 

a2ea2(£>_Zs) 

n 

u 

p2eO‘2[(D+L)-zs] 

a2eaA(D+L)~z’] 

(9) 


(10) 


(11) 


If  the  determinant  of  the  coefficient  matrix,  Det[C(k)],  vanishes,  then  the  solutions  of  the  Eq.  (9)  have 
poles  for  the  corresponding  horizontal  wavenumber  values.  The  poles  (pole)  correspond  to  additional 
modes  (mode)  in  the  wavenumber  spectrum. 
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3  Numerical  examples 

To  study  sound  propagation  in  the  Arctic  Sea  4  cases  axe  investigated.  The  measurements 
are  motivated  from  MIZEX  87  [4]  for  which  an  ice-layer  is  included.  It  is  necessary  to  em¬ 
phasize  that  the  roughness  and  other  inhomogeneous  properties  of  ice  are  neglected  in  these 
calculations.  For  each  case  some  comments  are  given  for  observations. 

Case  1:  3  layered  model,  Deep  waveguide  (vacuum-ice-fluid-fluid  waveguide) 

The  study  begins  with  an  ice- covered  sea  environment  for  which  a  half-space  fluid  is  in¬ 
cluded  to  neglect  the  effect  of  bottom  properties.  This  model  can  be  related  to  mathematical 
equations  presented  at  previous  section.  The  following  parameters  were  used  to  compute 
Green  s  function  and  normal  stress  shown  in  Fig.  1  and  2  for  two  different  receiver  depths: 

Table  1.1  :  Environmental  parameters  for  the  vacuum-solid-fluid-fluid  waveguide. 


ct(m/s) 

ca(m/s) 

p(gr/cm 3) 

depth  (m) 

Layer  0  (vacuum) 

oo 

Layer  1  (ice) 

ci  =3000 

D  =  4 

Layer  2  (water) 

c2  =1441 

II 

CO 

o 

Layer  3  (water) 

c3  =1457 

0 

00 

Where  c/  is  the  compressional  sound  velocity  and  cs  is  the  shear  sound  velocity  and  p  is 
density. 


Figure  1:  Results  for  Table  1.1  in  layer  2  (RD=18  m)  and  layer  3  (RD=122  m)  (a)  modulus  of  depth- 
separated  Green’s  function  and  (b)  transmission  loss  for  f=40  Hz,  N=4096,  contour  offset1  (e)=  0.053 
dB/A,  SD=7  m,  RD=18  m  (dashed  lines)  and  RD=122  m  (solid  lines). 


Figure  2.  Results  for  Table  1.1  (a)  modulus  of  depth-separated  Green’s  function  and  (b)  transmission 
loss  for  f=300  Hz,  N=16384,contour  offset  (e)=  0.013  dB/A,  SD=7  m,  RD=18  m  (dashed  lines)  and 
RD=122  m  (solid  lines). 


(1)  See  [[2],  p.35]. 
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Some  comments  are  given  for  observations: 

Green’s  function,  G2(k,z) 

-  In  Fig.  1  (a)  (/  =  40 Hz)  the  sharp  peak  is  a  discrete  pole  associated  with  the  shear 
sound  velocity  in  the  ice. 

-  In  Fig.  1  (a)  the  last  peak  ( k  ~  0.37  m~l)  is  an  additional  pole  which  appears  in  the 
evanescent  spectrum.  It  is  a  dramatic  shear  effect.  The  pole  corresponds  to  vanishing  of  the 
Det[C(k)\,  C  is  given  by  Eq.  (10).  It  belongs  to  the  class  of  leaky  waves.  The  peak  disappears 
when  the  receiver  is  moved  out  deeper  in  the  water  (RD=122  m),  because  in  Eq.  (7)  the  first 
exponential  term  becomes  very  small  as  the  receiver  depth  becomes  large. 

-  In  Fig.  2  (a)  (/  =  300 Hz)  The  sharp  peak  associated  with  sound  propagation  in  the 
ice  (shear  velocity)  is  a  discrete  pole.  The  peak  at  k  cz  0.6  m~l  is  a  part  of  the  continuous 
spectrum,  indicating  a  leaky  mode  associated  with  the  compressional  sound  velocity  in  the  ice. 

Normal  stress 

-  In  Fig.  1  (b)  for  z  =  18m  the  transmitted  field  in  the  layer  1  from  layer  2  can  be  both 
compressional  and  shear  fields  and  these  can  be  reflected  from  the  ice- vacuum  interface  and 
can  be  transmitted  again  as  compressional  field  in  the  layer  2.  The  interference  pattern  has 
been  dominated  by  the  evanescent  part. 

Case  2:  3  layered  model  (Vacuum,  Arctic  sound  velocity  profile) 


Table  1.2:  Environmental  parameters  for  the  Arctic  waveguide  1. 


c/i  [top]  (m/s) 

c/2  [bottom]  (m/s) 

p(gr/cm3) 

depth  (m) 

Layer  0  (vacuum) 

0 

0 

0 

oo 

Layer  1  (water) 

1439 

1441 

1.0273 

80 

Layer  2  (water) 

1441 

1456 

1.0278 

60 

Layer  3  (water) 

1456 

1461 

1.0279 

360 

where  c/i  is  the  compressional  sound  velocity  at  the  top  of  the  layer  and  c*2  compressional 
sound  velocity  at  the  bottom  of  the  layer  which  is  treated  by  SAFARI  as  a  fluid  with 
varying  linearly  with  depth.  D  is  the  ice  thickness  and  p  is  density. 

The  modulus  of  the  Green’s  function  and  normal  stress  for  the  environmental  data  given 
in  Tables  1.2  have  been  calculated  in  the  layer  1  and  2  by  SAFARI  and  shown  in  figures  3 
and  4  for  the  frequencies  40  and  300  Hz  respectively. 

As  one  can  observe  evanescent  field  is  important  only  for  low  frequencies  at  near  field  (in 
the  number  of  wavelength),  Fig.  3  (a). 

Waves  with  frequencies  lower  than  cutoff  frequency  f cutoff  =  - rj — “  can  no*  ProPagate 

along  the  layer  without  damping.  The  cutoff  frequency  in  layer  1  is  about  f cutoff  —  85  Hz 
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Figure  3:  Arctic  waveguide  1  (Table  1.2),  (a)  modulus  of  depth-separated  Green’s  function, 
(b)  transmission  loss  for  f=40  Hz,  N=4096,  contour  offset  (e) =0.053  dB/A,  SD=3  m,  RD=18 
m  (dashed  lines),  in  the  layer  1,  and  RD=122  m  (solid  lines),  in  layer  2. 


Figure  4:  Arctic  waveguide  1  (Table  1.2),  (a)  modulus  of  depth-separated  Green’s  function, 
(b)  transmission  loss  for  f=300  Hz ,  N=16384,  contour  offset  (e)=0.053  dB/A,  SD=3  m, 
RD=18  m  (dashed  lines),  in  the  layer  1,  and  RD=122  m  (solid  lines),  in  layer  2. 


and  in  layer  2  about  /cutoff  —  16  Hz.  Sound  propagating  with  frequencies  smaller  than 
cutoff  frequency  can  transmit  to  the  next  layer  and  if  the  source  frequency  is  larger  than  the 
cutoff  frequency,  then  the  sound  energy  can  be  channeled  in  that  layer.  In  Fig.  3  (b)  the 
sound  energy  leaks  from  layer  1  into  layer  2  and  is  channeled  and  concentrated  in  the  layer 
2.  Therefore  the  transmission  loss  at  the  receiver  depth  z  =122  m  is  less  than  that  of  in 
the  layer  1  (z  =  18  m).  This  observation  is  in  agreement  with  the  ambient  noise  recording 
measurements  made  in  Barents  Sea  during  SIZEX  92  [1], 

Case  3:  4  layered  model:  Vacuum,  Arctic  sound  velocity  profile  and  bottom 

The  modulus  of  Green’s  function  and  normal  stress  for  the  environmental  data  given  in 
Tables  1.3,  presented  in  case  3  (in  which  D,  ice  thickness,  is  zero),  have  been  calculated  and 
shown  in  figures  5,  6,  7  and  8  for  frequencies  40  and  300  Hz  respectively. 

In  the  next  part,  results  are  compared  with  case  4  for  which  an  ice-layer  is  included. 
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Case  4:  5  layered  model  (Vacuum,  ice,  Arctic  sound  velocity  profile  and  bottom) 


Table  1.3:  Environmental  parameters  for  the  Arctic  waveguide  2. 


on  [top]  (m/s) 

c;2[bottom](m/s) 

cs(m/s) 

p(gr/cm a) 

depth  (m) 

Layer  0  (vacuum) 

0 

0 

0 

0 

OO 

Layer  1  (ice) 

3000 

3000 

1600 

0.92 

D 

Layer  2  (water) 

1439 

1441 

0 

1.0273 

80 

Layer  3  (water) 

1441 

1456 

0 

1.0278 

60 

Layer  4  (water) 

1456 

1461 

0 

1.0279 

356 

Layer  5  (gravel*) 

1800 

1800 

0 

2.0 

OO 

*  A  compressions]  attenuation  (71  =  0.6  dB/A)  is  introduced  for  the  absorbing  bottom. 


The  modulus  of  Green’s  function  and  normal  stress  for  the  environmental  data  given  in 
Table  1.3,  with  an  ice  thickness  D  =  4m,  have  been  calculated  and  shown  in  figures  9,  10, 
11  and  12  for  frequencies  40  and  300  Hz  respectively.  Introduction  of  the  ice  has  reduced  the 
amplitude  of  the  sharp  peak.  There  are  no  significant  changes  in  the  corresponding  normal 
stresses  at  the  low  frequencies  comparing  with  the  last  model  in  which  there  was  no  ice  at 
the  top. 

As  frequency  increases  to  /  =  300  Hz,  normal  stresses  shown  in  figures  11  (b)  and  12  (b), 
have  been  damped  at  short  ranges,  comparing  with  figures  7  (b)  and  8  (b).  The  reason  is  that 
the  field  can  be  transmitted  in  the  ice  for  the  short  ranges.  The  transmitted  field  is  mostly 
absorbed  or  trapped  in  the  ice  layer. 

4  Summary  and  conclusions 

In  this  paper  SAFARI  is  applied  to  study  acoustic  propagation  in  an  ice-covered  ocean.  There 
are  several  attractive  features  of  SAFARI. 

•  The  code  is  numerically  stable  for  any  number  of  layers. 

•  The  field  at  a  number  of  receivers  can  be  determined  by  just  one  run  (not  shown  here)  . 

•  SAFARI  is  very  useful  for  the  treatment  of  interface  phenomena. 

However  some  numerical  limitations  have  been  observed. 

•  Numerical  difficulties  appear  when  the  frequency  becomes  large.  To  handle  this  the  number 
of  sampling  points  must  be  made  large, 

•  Some  very  short-range  effects  are  neglected  in  SAFARI,  because  of  applying  the  asymptotic 
approximation  of  the  Hankel  functions  for  the  wavenumber  integrand.  This  results  in  oscilla¬ 
tions  in  the  normal  stress  curves  at  very  short  ranges. 

The  study  is  performed  in  a  progression  from  simple  environments  to  more  complicated  mod¬ 
els  that  simulate  ocean  Arctic  environments.  Each  model  introduces  greater  complexity  so 
that  the  importance  of  the  different  parameters  can  be  isolated.  The  analytic  integral  expres¬ 
sions  for  the  field  for  a  simple  3  layered  model  are  obtained  by  using  boundary  conditions.  In 
these  calculations  only  compressional  monochromatic  point  sources  have  been  used.  In  each 
experiment  both  the  source  and  receivers  are  localized  in  the  water  below  the  ice  interface. 
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The  computations  by  SAFARI  are  performed  for  source  frequencies  at  40  and  300  Hz.  The 
main  observations  and  results  from  the  numerical  study  are: 


•  As  the  separation  between  the  receiver  and  the  interface  increases,  either  by  increasing 
of  the  frequency  or  for  a  fixed  frequency  by  moving  the  receiver  down  away  from  the  inter¬ 
face,  the  exponential  decay  makes  the  evanescent  portion  diminished. 

•  Sound  will  be  transmitted  into  an  ice  layer  as  compressional  and  shear  waves.  There  is 
a  competition  between  compressional  and  shear  waves  at  the  ice-water  boundary. 

•  For  the  incident  angles  less  than  the  critical  angle  (corresponding  to  the  compressional 
waves),  there  are  only  compressional  waves  in  the  ice  layer.  As  the  range  increases  the  shear 
waves  transmission  coefficient  increases  while  the  transmission  coefficient  for  the  compres¬ 
sional  waves  decreases. 

•  At  a  particular  range  the  transmission  coefficients  for  the  compressional  waves  and  shear 
waves  become  equal  and  interfere  destructively.  This  phenomena  can  be  observed  as  the 
regularly  spaced  frequency  resonances  in  reflection  loss  patterns  and  as  a  local  minimum  in 
transmission  loss  patterns. 


•  The  subsequent  shear  waves  are  trapped  in  the  ice  and  even  a  small  shear  attenuation 
eventually  removes  all  the  energy  in  the  shear  field.  The  trapped  compressional  waves  in  the 
ice  layer  are  only  important  at  short  ranges. 

•  The  elastic  properties  of  the  ice  have  a  significant  effect  on  the  reflection,  although  these 
properties  do  not  affect  the  frequency  separation  of  the  peaks  in  the  reflection  loss  pattern. 
The  strengths  of  the  mode  conversions  at  the  water-ice  boundary  depend  on  the  ice  compres¬ 
sional  and  shear  sound  speed. 
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Figure  5:  Results  for  Arctic  waveguide  2  (Table  1.3),  for  D— 0,  (a)  modulus  of  depth-separated 
Green’s  function  in  layer  2,  |  G2{k,z )  |,  and  (b)  Transmission  loss  for  f=40  Hz,  N=4096,  contour 
offset(e)=0.053  dB/A,  SD=3  m,  and  RD=18  m. 


Figure  6:  Results  for  Arctic  waveguide  2  (Table  1.3),  for  £>=0,  (a)  modulus  of  depth-separated 
Green’s  function  in  layer  3,  |  Gz{k,  z)  |,  (b)  transmission  loss  in  layer  3  for  f=40  Hz,  N=4096, 
contour  offset  (e)=0.053  dB/A,  SD=3  m,  RD=122  m  . 


Figure  7:  Results  for  Arctic  waveguide  2  (Table  1.3),  for  D- 0,  (a)  modulus  of  depth-separated 
Green’s  function  in  layer  2,  |  G2(k,z)  j,  (b)  transmission  loss  in  layer  2  for  f=300  Hz, 
N=16384,  contour  offset  (e)=0.013  dB/A,  SD=3  m,  RD=18  m. 


Figure  8:  Results  for  Arctic  waveguide  2  (Table  1.3),  for  D= 0,  (a)  modulus  of  depth-separated 
Green’s  function  in  layer  3,  |  G3(k,z )  |,  (b)  Transmission  loss  in  layer  3  for  f=300  Hz, 
N=16384,  contour  offset  (e)=0.013  dB/A,  SD=3  m  and  RD=122  m. 
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Figure  9:  Results  for  Arctic  waveguide  2  (Table  1.3),  for  D= 4  m  (ice  thickness),  (a)  modulus 
of  depth-separa  ted  Green’s  function  ,  |  G2(k,  z)  |,  (b)  transmission  loss  for  f=40  Hz,  N=4096, 
contour  offset (e) =0.053  dB/A,  SD=7  m,  and  RD=18  m. 


Figure  10:  Results  for  Arctic  waveguide  2  (Table  1.3),  for  D= 4  m  (ice  thickness),  (a)  modulus 
of  depth-separated  Green’s  function  in  layer  3,  |  G3(k,z)  |,  (b)  transmission  loss  in  layer  3  for 
f=40  Hz,  N=4096,  contour  offset  (e)=0.053  dB/A,  SD=7  m,  RD=122  m. 


Figure  11:  Results  for  Arctic  waveguide  2  (Table  1.3),  for  D= 4  m  (ice  thickness),  (a)  modulus  of 
depth-separated  Green’s  function  in  layer  2,  |  G2{k,z )  |,  (b)  transmission  loss  in  layer  2  for  f=300  Hz, 
N=16384,  contour  offset  (e)=0.013  dB/A,  SD=7  m,  RD=18  m. 


Figure  12:  Results  for  Arctic  waveguide  2  (Table  1.3),  for  D= 4  m  (ice  thickness),  (a)  modulus  of 
depth-separated  Green’s  function  in  layer  3,  |  G3(k,z )  |,  (b)  transmission  loss  in  layer  3  for  f=300  Hz, 
N=16384,  contour  offset  (e)=0.013  dB/A,  SD=7  m  and  RD=122  m. 
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Abstract 

This  numerical  study  describes  reflection  of  plane  acoustic  waves  from  an  ice-layer  overlying 
a  water  halfspace.  Compressional  reflection  coefficient  and  reflection  loss  for  a  simple  two 
layered  model  are  investigated.  Results  illustrated  the  dependence  of  reflection  coefficient 
on  ice  thickness  and  source  frequency. 

1  Introduction 

The  discontinuity  in  the  interface  between  two  materials  produces  a  significant  influence  on 
the  wave  propagating  through  the  medium.  The  reflection  coefficient  is  a  mathematical  formu¬ 
lation  of  the  physical  effect  of  the  boundary  on  the  incident  wave  field  and  can  be  mentioned 
as  a  filter  on  the  wavenumber.  A  large  amount  of  literature  deals  with  acoustic  reflection 
from  the  bottom  or  sea  floor.  The  sound  reflection  from  the  ice-covered  sea  environment 
can  be  treated  similarly.  Ice  can  be  considered  as  an  elastic  medium,  which  means  an  inci¬ 
dent  compressional  plane  wave  can  be  reflected/transmitted  as  both  compressional  and  shear 
waves.  The  reflection  coefficient,  the  ratio  of  the  amplitude  of  the  reflected  plane  wave  to  the 
amplitude  of  a  plane  wave  incident  on  the  ice  overlying  the  sea,  is  an  important  factor  for 
investigation  of  the  effect  of  the  ice-cover  on  acoustic  propagation.  The  roughness  properties 
of  ice  are  neglected  in  this  study. 

This  paper  focuses  on  the  ratio  of  the  amplitude  of  the  direct  reflected  compressional  plane 
wave  to  the  amplitude  of  an  incident  compressional  plane  wave,  as  shown  in  Fig.  1: 


Z 

Figure  1:  Reflection  of  an  incident  compressional  (P)  plane  wave  from  interfaces  of  a  vacuum-ice-water 
environment,  PP:  direct  reflected  P  wave. 


Hydrological  measurements  in  the  Marginal  Ice  Zone  (MIZ)  of  the  Greenland  Sea  reveals  that 
the  sound  speed  profile  for  the  ice-covered  Arctic  Sea  can  be  considered  range-independent. 
SAFARI,  which  provides  accurate  transmission  loss  values  for  propagation  in  range-independent 
multi-layered  horizontally  stratified  environment,  can  be  a  useful  tool  to  study  the  reflection 
coefficient  for  the  ice-covered  region  of  the  Arctic  Sea. 

In  order  to  provide  a  basis  for  an  initial  investigation,  the  reflection  coefficient  formula,  for 
direct  reflected  P  wave  for  a  simple  2  layered  environment  (vacuum-solid-fluid  halfspace)  is 
presented  in  Sec.  2  .  A  parameter  study  of  effects  of  compressional  and  shear  properties  of 
ice  on  the  reflection  loss  is  given  by  a  numerical  example  in  Sec.  3.  Dependency  of  reflection 
coefficient  on  incident  angle,  frequency  and  ice  thickness  is  studied  in  Sec.  3.2.  Finally  the 
results  are  summerized  in  Sec.  4.  The  simulations  were  performed  by  SAFARI  [1], 


2  Mathematical  formulation 


The  compressional  reflection  coefficient,  RpP,  for  the  3  layered  environment  shown  in  Fig.  1 
(vacumm-solid-fluid)  can  be  expressed  as: 


Rpp  =  e  l<f>,  (j>  =  2 arctan[ 


(MZi)Zs 


(MZi)2  -  (NZi)2i 


where 


M  =  'z[\-Z^cos^s)2cotg{p)  +  Zssin(26s)2cotg(q)\, 

N  _  Z2Cos(26s)2  +  Zssin(26s)2 

Z\sin(p )  Z\sin{q)  ' 


p  =  k2dcos(02),  q  =  ksdcos(6s), 


_  P 1C*  Z  _  P2C2  7  P2CS  _  piCq 

cos(6\Y  2  cos(02)’  s  cos{Os)'  3  cos{03)' 


(1) 


(2) 

(3) 

(4) 


where  pup2  and  p3  are  densities  of  the  layers  1,2  and  3  (px=  0  for  the  vacuum),  Zx  and  Z3 
impedances  for  the  layer  1  (  Z\  —  0  for  vacuum)  and  3,  Z2  and  Zs  are  compressional  and  sheax 
impedances  for  layer  2,  6\  is  the  incident  angle,  03  is  the  transmitted  angle  in  the  layer  3,  O2 
and  6S  are  the  transmitted  angles  in  the  layer  2  corresponding  to  the  compressional  and  shear 
fields,  k2  and  ks  are  the  wavenumbers  in  layer  2  (ice)  corresponding  to  the  compressional  and 
shear  fields  [2]  .  The  reflection  coefficient  is  an  oscillatory  function  of  p  and  q. 


The  complex  reflection  coefficient  Rpp  can  be  expressed  by  modulus  |  RpP  \  and  phase 

*an  1  JRe[RpP}\ '  underwater  acoustics,  the  modulus  is  usually  represented  by  the  corre¬ 
sponding  reflection  loss: 


Rloss  —  20 log  |  RpP  |  (5) 

In  the  next  section  a  numerical  study  will  be  performed  on  the  compressional  reflection 
coefficient,  RpP ,  and  reflection  loss  for  a  simple  vacuum- ice- water  environment  in  which  a 
monochromatic  point  source  localized  just  below  the  ice  interface  in  the  water. 
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3  Numerical  example:  ice-covered  deep  waveguide 

A  simple  2  layered  environment  was  used  to  study  an  ice-covered  region  of  Arctic  Sea.  Even 
if  the  chosen  simulated  environment  is  very  simple,  remarkable  results  are  observed. 


Table  1.1:  Environmental  data  for  vacuum-ice-water  model. 


ci  {m/s) 

cs  (m/s) 

7 1  (dB/A) 

7s  (dB/A) 

p(gr/cmA) 

depth  (m) 

Vacuum 

0 

0 

0 

0 

0 

oo 

Ice 

3000 

1600 

1. 

2.5 

0.92 

d=  1,  2,  4 

Water 

1500 

0 

0 

0 

1.0 

00 

Where  cj  is  the  compressional  sound  velocity  and  cs  is  the  shear  sound  velocity,  p  is  den¬ 
sity  and  ji  and  7 s  are  compressional  and  shear  attenuation  coefficients  for  ice  respectively. 

Fig.  2  presents  the  modulus  of  the  reflection  coefficient  versus  incident  angle  which  con¬ 
firms  the  dependency  of  reflection  coefficient  on  frequency,  /,  and  ice  thickness,  d. 

Frequency =  5  KHz  Frequency =  20  KHz 


Figure  2:  Modulus  of  the  compressional  reflection  coefficient  for  Table  1.1  for  ice  thickness  4  m  (solid 
curve),  2  m  (dashed)  and  1  m  (dotted). 


3.1  Interesting  cases  for  vacuum- ice- water  model 

In  this  part  reflection  loss,  given  by  Eq.  (5),  is  studied  for  some  important  angles  or  different 
sound  velocities.  Reflection  loss  for  the  following  angles  is  calculated  and  is  shown  in  Fig.  3: 

1)  0  =  0°,  normal  incident. 

2)  0  =  30°  (the  critical  angle  due  to  the  compressional  sound  velocity  in  the  ice). 

3)  0=41.5°  ,  for  which  the  pure  shear  wave  transmission  occurs. 

4)  0  =69.5°  (the  critical  angle  due  to  the  shear  sound  velocity  in  the  ice). 


23 


Incident  Angle*  0.0  (deg.) 


Incident  Angle*30.0  (deg) 


Frequency{Hz) 

0  1000  2000 


0  1000  2000 

Frequency  (Hz) 


Incident  Angle* 41  £  (deg) 


Frequency(Hz) 


Frequency(Hz) 


/nc<d«nt  Angles  6  9. 6  (deg) 


Frequency  (Hz)  Frequency(Hz) 


Figure  3.  Reflection  loss  ( Rioss  —  20 log  \  Rpp  |)  for  vacuum-ice-water  model,  ice  thickness 
d= 4  m  (solid  curve),  2  m  (dashed)  and  1  m  (dotted). 


The  frequency  of  the  peaks  corresponds  to  2dcos(es )  >  so  when  the  ice  thickness  increases  the 
spacing  of  the  peaks  decreases.  These  peaks  appear  when  the  reflection  coefficient  becomes 
very  small,  and  the  transmitted  waves  are  trapped  in  the  ice  layer. 

To  illustrate  how  the  variation  of  different  parameters  that  describe  the  physical  proper¬ 
ties  of  environment  can  affect  the  reflection  loss,  the  following  cases  are  considered 

a)  Variation  of  the  compressional  velocity  in  the  ice  is  studied  in  Fig.  4  (a)  (co=  3500 
3200  and  3000  (m/s)). 

b)  Variation  of  the  shear  sound  velocity  in  the  ice  is  studied  in  Fig.  4  (b)  (c5=1800,  1700  and 
1600  (m/s)). 

The  reflection  loss  is  sensitive  to  the  compressional  sound  velocity  in  the  ice  for  6  <  9^  =  30°, 
where  compressional  waves  dominate.  For  6  >  6cr\ ,  where  the  shear  waves  dominate,  reflec¬ 
tion  loss  increases  by  increasing  of  the  shear  sound  velocity  in  the  ice. 
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(a)  cs=1600  m/s,  variation  of  c2:  3500  (solid  (b)  c2=3000  m/s,  variation  of  cs:1800  (solid 

lines),  3200  (dashed  lines)  and  3000  m/s  (dot-  lines),  1700  (dashed  lines)  and  1600  m/s  (dot¬ 
ted).  ted). 

Figure  4:  Reflection  loss  ( RdB  =  -20 log  \  Rpp  |)  for  Table  1.1  (vacuum-ice-water) ,  Fre- 
quency=10  KHz  and  ice  thickness=4  m. 

3.2  Numerical  results 

For  a  vacuum-ice- water  environment,  as  shown  in  Fig.  1,  the  reflection  coefficient  is  a  function 
of  the  incident  angle,  source  frequency  and  the  ice  thickness  (R —R(9,f,d)).  Contour  plots 
can  be  used  to  demonstrate  the  effect  of  frequency  for  different  incident  angles  for  a  fixed  ice 
thickness  on  reflection  coefficient.  Contouring  is  performed  for  the  samplings  A 9  —  0.5°  and 
A /  =  200 .  Fig.  5.  demonstrates  contour  plots  for  the  modulus  of  the  reflection  coefficient  for 
the  parameters  given  in  Table  1.1. 

The  contour  plots  can  be  discussed  by  the  multiplication  of  the  source  frequency  and  ice 
thickness  term  ( df )  or  by  a  relation  between  d  and  A2,  where  A2  is  the  wavelength  of  the 
compressional  waves  in  the  ice  (m). 

-  For  df  <  200  (d  <  y|),  total  reflection  occurs  for  all  incident  angles. 

-  For  200  <  df  <  500  (d  <  ^),  there  is  a  phase  shift  at  the  critical  angle  (9cr\  =  30°). 
No  shear  waves  have  been  exited  in  ice  yet. 

-  For  <  df  <  800  (d  <  transmitted  shear  waves  can  be  observed  for  9  >  9cr\  —  30°. 

-  For  df  >  800  (d  >  5^5): 

There  are  oscillations  for  9  <  9cr  1  =  30°,  where  the  compressional  waves  dominate  and 
oscillations  reduces  as  df  becomes  larger. 
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The  reflection  coefficient  becomes  smaller  and  smaller  for  41°  <  9  <  55°,  where  shear  waves 
dominate  and  corresponds  to  the  high  reflection  loss  region.  In  this  interval  interface  waves 
of  the  Scholte  or  Stonely  [2]  type  are  excited. 

-  For  df  >  4000  ( d  <  ^l2-)  and  9  >  6cr 2  ~  70°,  reflection  coefficient  does  not  vary  any 
more. 


-  For  df  >  12000  (d  >  4A2),  the  refection  coefficient  in  the  interval  41°  <  6  <  55°  has 
become  very  small  and  stabilized. 


-  The  reflection  coefficient  becomes  very  small  when  the  frequency,  the  ice  thickness  and 

the  incident  angle  confirm  to  the  relation  f  = - .  CjC] 

J  2  d,Jc\-clsin(e) 

4  Summary 


Compressional  waves  are  divided  to  the  shear  waves  and  compressional  waves  at  the  boundary 
of  the  water-ice  interface.  These  shear  waves  can  be  reflected  from  the  ice- vacuum  boundary 
and  converted  back  to  the  compressional  waves  into  water  through  the  ice-water  boundary. 
This  conversion  can  be  observed  especially  when  the  frequencies  and  the  angles  and  the  ice 
thickness  coincide  with  shear  wave  resonances  in  the  ice  layer  (/  =  -2dcffe  y),  which  causes 

the  noticeable  peaks  in  the  reflection  loss  (see  Fig.  3).  The  mode  conversion  coefficients  are 
sensitive  to  the  compressional  and  the  shear  speed  in  the  ice. 


A  study  is  performed  for  different  incident  angles  and  frequencies  and  ice  layer  of  thick¬ 
ness  1,  2  and  4  m  ,  d=  |,  A,  2A,  which  can  be  an  interesting  area  for  some  further  experiments. 

At  low  frequencies  and  low  angles,  the  transmitted  compressional  wave  is  evanescent,  only 
the  direct  reflection  is  present  and  the  field  decays  exponentially  in  the  ice  separating  from 
the  boundary.  If  the  ice  is  very  thin,  these  compressional  waves  can  be  excited  as  surface 
waves. 


The  shear  properties  of  the  ice  dominate  when  the  angle  of  incidence  exceeds  the  critical  angle 
related  to  the  compressional  sound  velocity  in  ice  (approximately  6  =  30°)  and  the  reflection 
loss  is  very  high  for  the  interval  30°  <  9  <  69.5°.  In  this  case  sound  can  penetrate  into  the  ice. 

The  reflection  coefficient  for  the  incidence  angles  larger  than  the  critical  angle  regarding 
to  the  shear  sound  velocity  in  ice  (approximately  9  =  70°)  is  very  complicated  since  both 
transmitted  angles  are  complex. 

In  the  vacuum-ice-water  environment,  ice  thickness  is  not  important  for  high  frequencies 
and  large  incident  angles. 
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SUMMARY 

The  purpose  of  this  study  has  been  to  develop  a  numerical  model  for  sound  propagation  in  both  air  and 
ground  over  long  distances.  Wave  propagation  in  both  air  and  ground  is  modeled  by  a  set  of  partial 
differential  equations  for  wave  propagation  in  a  porous  medium  with  a  rigid  frame.  By  introducing 
Chebyshev  spectral  collocation  in  the  spatial  variable  the  equations  are  transformed  to  a  set  of  ordinary 
differential  equations  in  time.  Domain  decomposition  is  used  to  divide  the  computational  domain  into 
subdomains  of  manageable  sizes.  All  subdomains  are  homogenous  in  the  sense  that  they  contain  only 
ground  or  air.  Some  numerical  examples  are  calculated  in  two  dimensions  and  the  results  are  compared 
to  analytical  results  and  to  some  results  from  outdoor  measurements.  These  comparisons  show  the  low 
dispersion  and  dissipation  of  spectral  collocation,  and  that  the  suggested  model  gives  good  results  for 
sound  propagation  over  long  distances. 

THE  DIFFERENTIAL  EQUATIONS 

Sound  propagation  in  a  porous  medium  with  a  rigid  frame  is  modeled  by  the  linear  Eqs.  (1.1)  and  (1.2). 


Vp  +  p' —  +  Rv  -  0 
dt 

— +  p'c'2V-v  «0 
dt 


(1.1) 

(1.2) 


Here  p  is  the  air  pressure,  v  a  vector  containing  the  particle  velocity  components  and  R  the  flow 
resistivity.  Furthermore  is  p'  the  equivalent  density  and  c  the  equivalent  sound  speed  for  the  porous 
medium  and  they  are  given  by  Eq.  (1.3)  and  (1.4)  respectively. 


(1.3) 


(1.4) 


Here  is  p<>  the  air  density,  ks  the  structure  factor,  cp  the  porosity  and  c  the  sound  speed  in  the  air. 
SPECTRAL  COLLOCATION 

If  a  two-dimensional  Cartesian  coordinate  system  is  chosen,  Eq.  (1.1)  and  (1.2)  can  be  written  as  Eq. 

(1.5). 


3U  dU  4  dU 

—  +  A,  —  +  A,  —  +  CU  =  0 
dt  dx  "  dy 


The  vector  U  is  defined  by  U  =  [u  v p]  with  u  and  v  being  the  horizontal  and  vertical  particle  velocity 
components  and p  the  air  pressure.  Introducing  Chebyshev  spectral  collocation  in  the  spatial  variable 
leads  to  a  discretization  in  space  and  thereby  a  transformation  of  our  PDEs  to  ODEs  in  time. 


3  1  R  n 

— u  +  — Dxp  +  — u  =0 

dt  p'  p' 

d  1  «  R  n 

—  v  +  —  D  p  +  -v  =  0 
3t  P  P 

— p  +  P'c'2(d  u  +  D  v)  =  0 
dt 


(1.6) 

(1.7) 

(1.8) 
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The  matrixes  u,  v  and  p  contain  the  velocity  components  and  pressure  values  at  the  discretization 
points,  which  are  called  collocation  points  in  the  spectral  collocation  method.  Matrixes  Dx  and  Dy  are 
the  derivative  matrixes  in  the  two  spatial  directions  (see  Canuto,  Hussaini,  Quarteroni  and  Zang  1988 
for  more  details  about  the  derivative  matrixes).  A  second  order  explicit  Runge-Kutta  method  was  used 
in  the  time  integration  of  the  ODEs. 

DOMAIN  DECOMPOSITION 

Domain  decomposition  was  introduced  to  be  able  to  calculate  sound  propagation  for  large 
computational  domains.  The  computational  domains  were  divided  into  a  set  of  smaller  non-overlapping 
subdomains.  Computational  domains  that  contained  both  air  and  ground  were  divided  into  subdomains 
containing  only  air  or  ground,  so  in  a  sense  the  subdomains  were  homogenous.  Satisfaction  of  the 
differential  equations  had  to  be  ensured  on  the  boundaries  between  subdomains  too.  Therefore  a 
correctional  method,  where  the  boundary  values  were  corrected  between  each  timestep  in  the  Runge- 
Kutta  time  integration,  was  introduced.  The  corrected  values  were  calculated  oh  the  basis  of  the 
physical  boundary  conditions,  continuous  normal  particle  velocity  and  air  pressure  ,  and  implemented 
via  characteristic  boundary  conditions  (Bj0rhus  1995  and  2000). 

OPEN  BOUNDARIES 

Reflections  of  the  propagating  signal  from  the  outer  boundaries  of  the  computational  domain  have  a 
tendency  to  contaminate  the  solution.  Such  reflections  should  therefore  be  kept  as  small  as  possible.  An 
approximation  to  the  exact  open  boundary  conditions  was  therefore  implemented  via  characteristic 
boundary  conditions  at  the  outer  boundaries.  The  chosen  approximation  was  quite  good  for  normal 
incidence,  but  unfortunately  it  didn’t  perform  nearly  as  good  when  the  angle  of  incidence  increased. 

MAPPING  OF  DOMAINS 

Because  spectral  collocation  requires  square  domains,  some  special  attention  had  to  be  made  when 
domains  of  more  complex  shapes  were  used.  This  problem  was  solved  by  performing  a  mapping  of 
complexly  shaped  domains  onto  a  square.  The  mapping  lead  to  a  transform  of  the  spatial  variables. 

NUMERICAL  EXPERIMENTS 

Some  numerical  experiments  were  calculated  to  test  the  proposed  model.  The  calculations  were  done  in 
two  dimensions.  A  computational  domain  with  66  subdomains  (3x22)  and  the  dimensions  15x110 
meters  (height  x  length)  was  used.  There  were  51  collocation  points  in  each  direction  in  each 
subdomain.  This  gave  a  resolution  of  3.44  points  per  wavelength  for  a  1000  Hz  signal.  A  point  source 
was  implemented  by  trigging  one  collocation  point  with  an  air  pressure  that  varied  with  time.  This 
point  source  was  set  to  radiate  a  bandlimited  signal  between  approximately  100  and  900  Hz.  The  power 
spectral  density  of  the  source  signal  is  shown  in  Fig.  1. 

The  sound  speed  and  damping  of  the  signal  in  free  field  was  calculated  from  a  simulation  with  only  air 
domains.  The  sound  speed  between  5  and  80  meters  from  the  source  is  shown  in  Fig.  1.  The  damping 
of  the  signal  over  the  same  distance  is  shown  in  Fig.  2. 


Fig.  1 .  Sound  speed  between  5  and  80  meters  from  the  source.  Fig.  2.  Damping  between  5  and  80  meters  from  the  source. 

In  these  experiments  a  sound  speed  of  344  m/s  was  chosen  for  the  air.  Fig.  1  shows  that  the  sound 
speed  in  the  calculations  is  highly  accurate  within  the  frequency  band  of  the  source  signal.  In  other 
words  there  is  hardly  any  dispersion  at  distances  up  to  at  least  80  meters.  For  two  dimensions  a 
damping  of  12  dB  is  expected  between  5  and  80  meters  from  the  source. 
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Fig.  2  shows  that  the  damping  in  the  numerical  experiments  is  in  good  accordance  with  the  expected 
damping.  There  is  however  a  dip  centered  at  900  Hz.  This  dip  is  most  likely  caused  by  interference 
between  the  direct  wave  and  a  reflection  from  the  boundaries  of  the  computational  domain.  This  shows 
the  importance  of  good  open  boundary  conditions  at  the  outer  boundaries  of  the  computational  domain. 

An  approximation  to  the  ground  profile  shown  in  Fig.  3  was  implemented  based  on  measurements  done 
at  Gran&sen  in  Trondheim  (Storeheier  1999). 


Ground  profile 


Fig.  3.  Measured  ground  profile, 

The  parameters  in  the  ground  were  chosen  based  on  measurements  done  at  this  site.  At  zero  meters  a 
point  source  was  placed  1.5  meters  above  the  ground  (0,4.5).  The  source  signal  was  the  same  as  for  the 
experiment  described  above.  The  resulting  air  pressure  was  registered  at  a  point  2  meters  above  the 
ground  and  100  meters  from  the  source  (100,6.5).  The  air  pressure  relative  to  the  expected  free  field 
value  for  two  dimensions  (insertion  loss)  was  calculated  and  compared  to  measurements  done  at  the 
site. 


Fig.  4.  Air  pressure  relative  to  free  field  at  100  meters  from  the  source. 

The  results  from  the  calculations  with  spectral  collocation  follow  the  measurements  quite  nicely  except 
for  the  fact  that  the  destructive  interference  at  250  Hz  isn ’t  large  enough. 

CONCLUSION 

Spectral  collocation  is  well  known  for  its  low  dispersion  and  dissipation.  The  results  above  have 
confirmed  this  and  shown  that  spectral  collocation  is  a  useful  method  for  calculating  sound  propagation 
over  long  distances.  The  only  major  advantage  is  the  fact  that  the  calculations  are  very  time  consuming. 
It  took  40  hours  to  calculate  the  results  shown  in  Fig.  4.  The  calculations  were  done  on  a  PC  with  a  733 
MHz  Pentium  III  processor  and  512  Mb  of  internal  memory.  There  are  however  ways  to  reduce  the 
computational  time  and  this  should  be  investigated  further. 
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Introduction 

The  traditional  split-beam  echo  sounder  measures  the  phase  angle 
between  the  signals  received  on  two  transducer  halves.  With  four 
transducer  sections,  arranged  as  four  quadrants,  the  direction  to  the 
target  can  be  measured  both  longitudinal  and  transverse.  This  has  been 
the  principle  in  all  Simrad  split-beam  transducers  for  target  strength 
measurements;  from  the  first  ES380  in  1983  to  the  latest  200  kHz 
composite  transducer  ES200-7  introduced  last  year. 

However,  three  sections  would  be  enough  to  provide  the  necessary 
information;  and  a  new  split-beam  transducer  ES30-10  with  three 
sections  is  now  introduced.  The  frequency  is  38  kHz,  and  the  beamwidth 
is  10  degrees.  The  three  sections  are  identical,  with  120  degrees 
rotation  from  one  to  another. 


ES38-10 

with  three  sections 


split-beam  transducer 
with  four  sections 

bird’s-eye  view 


Angle  measurement 

For  a  traditional  split-beam  transducer  configured  with  two  halves,  the 
angle  u,  between  the  acoustic  axis  and  the  direction  to  the  target  is 
derived  from 

<p  =  kd  sin  u  (1) 

where  <p=phase  angle  between  the  signals  received 
on  the  two  transducer  halves 
/c=wave  number=2jr/A 
A=wavelength 

d=distance  between  the  centres  of  the  two  transducer  halves 
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Normally  the  angle  u  Is  small,  so  that  sin u  can  be  approximated  by  u, 
and  we  have 

(p-  kd  u 

The  product  A-kd  is  called  the  angle  sensitivity,  and  the  relation 
between  the  phase  angle  and  the  physical  angle  to  the  target  can  be 
written  as 


q>-Au  (2) 

In  Simrad  split-beam  echo  sounders  the  phase  angle  cp  is  detected  in  the 
receiver  and  transferred  to  the  main  processor,  where  the  angle  u  is 
calculated  as  u-q>/A 


For  a  transducer  with  three  sections  the  mathematics  is  a  little  more 
complicated.  The  acoustic  centres  of  the  three  sections  have  the  co¬ 
ordinates: 


~V  0) 

«2=(Y«  0,  (3) 

e3=(0  a  0) 

in  a  right  handed  xyz-system  with  the  y-axis  pointing  ahead  on  the 
vessel  and  the  z-axis  pointing  vertical  down,  a  is  the  distance  from 
origo  to  one  of  the  acoustic  centres.  For  the  new  38  kHz  transducer  a  is 
66  mm. 
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s  is  a  unit  vector  in  the  direction 
towards  the  target 

s  =  (cos  vx  COSVj,  cosvx) 

where  vx  vy  and  vt  are  the  angles 

from  s  to  the  axes  of  the  co-ordinate 
system. 

The  echo  signals  received  on  the 
three  transducer  sections  have 
different  phase  angles  caused  by  the 
different  travel  distances  to  the  three 
acoustic  centres.  With  origo  as  the 
reference  point,  the  extra  travel 
distance  to  a  point  e  is  the  inner 
product  e  s.  The  phase  angles  at  the 
acoustic  centres,  relative  to  origo,  are 


q>io  =  k  &i"S 

<p2o  -  k  e2s  (4) 

(p30  ~  k  Bj-S 


The  three  phase  angles  are  not  independent.  The  relation  between  them 
is 


<pio+<P2o+<P3o  -  k(ei+  e2+  e3)  s  =  0  (5) 

The  angles  vx  and  vy  can  be  derived  from  the  equations  (4)  and  (3) 


,  i  /  V3  1  , 

<p10=k  et  s=/c( - acosvx  acosvy) 

2  2 

/  /  i  , 

(p2o-k  e2s=k( — a cosvx  — arcosv^) 

2  2 


(6) 


(p30~k  e3  s=kacosvy 


which  lead  to 


cosv'=S>^°> 

C0S Vy~^ca  ^2<Pvi  ~ ^10  “ ^20 ) 


(7) 


The  echo  sounder  screen  presents  a  horizontal  cross-section  of  the 
beam,  with  fish  echoes  marked  as  dots  at  their  respective  positions  in 
the  beam.  The  unit  vector  s  is  pointing  in  the  direction  to  the  target,  and 
the  projection  of  s  onto  the  xy-plane  is  used  for  the  screen  presentation. 
The  projection  of  s  has  the  co-ordinates  (cosvx  cosv^,).  These  cosines 

have  to  be  calculated  for  every  depth  sample. 
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The  signals  from  the  three  transducer  sections  goes  to  three  receiver 
channels  and  comes  out  as  three  digital  complex  numbers,  G*  ,  G2  and 
G3.  In  the  equations  (7)  there  are  three  phase  angles.  They  can  be  found 
by  an  arctan(lm/Re).  However,  the  arctan  is  a  time  consuming  function, 
and  it  may  be  worthwhile  to  reduce  the  number  of  phase  angles  from 
three  to  two  with  the  following  manipulation.  We  define 


^31=^30-^10 

(piZ-<p30-(p20 

and  the  equations  (7)  can  be  written  as 

C°SV*  =~&ka^1  ~<P^ 
cosvy  =—(^3i  +<Pu) 


(8) 


(9) 


A  phase  angle  difference,  as  in  equation  (8)  can  be  deduced  from  the 
phase  angle  in  the  product  of  the  first  complex  number  and  the  other 
conjugated: 


(psi  -  Arg(G3  Gr*)+n-27t  (10) 

<p32  =  Arg(G3-G2*)+/7-27c 

where  *  means  complex  conjugate,  and  Arg  (the  argument)  is  a  phase 
angle  between  -it  and  it. 

For  fish  in  the  main  lobe  of  this  transducer  and  between  the  -10  dB 
points,  the  phase  angle  difference  in  equation  (8)  is  less  than  160 
degrees,  which  means  that  there  is  no  2%  ambiguity;  n=0. 

In  existing  Simrad  split-beam  echo  sounders,  the  complementary  angle 
to  v  is  often  used 

ux=7r/2-vx  sinMI=cosvI 

uy  =7r/2-vy  sin  uy  =co  sVy 

ux  and  Uy  are  the  angles  from  s  to  respectively  the  yz-plane  and  the  xz- 

plane.  Most  split-beam  transducers  are  narrow  beam  transducers,  so 
that  the  approximation  sinu=u  is  allowable  in  the  main  beam;  and  we 
have 


Ux  —  (^31  ^32 )  (11) 

u,  =^-(9,1+9^ 
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The  angle  sensitivity  A  is  much  used  in  existing  split-beam  echo 
sounders  related  to  transducers  with  four  sections.  For  a  transducer  with 
three  sections  the  use  of  angle  sensitivity  is  more  dubious.  However,  in 
order  to  make  use  of  existing  programs  in  the  echo  sounder  processor, 
we  state  a  similar  definition  of  the  angle  sensitivity 


A-kd 


where  d  is  the  distance  between  the  centres  of  two  of  the  sections 
c/=  V3 -  a 

The  signal  processor  in  the  receiver  calculates  two  ‘phase  angles’,  <px 
and  <py,  afthough  at  least  one  of  them,  <py  is  fictive  regarded  as  a  phase 
angle.  With  the  aid  of  equation  (1 1)  we  have 

<Px  =Aux=kdux=(fhr  (fh2  (12) 

(Py  ^JI  +  ^32)/  V3 


<px  and  <py  are  transferred  to  the  main  processor,  where  the  angle  to  the 
target  is  calculated  as 

ux-q>x/A  (13) 

uy=<p/A 

For  this  transducer  vl=V3far=18.3 

The  coordinate  system  used  here  has  the  x-axis  pointing  port  and  the  y- 
axis  pointing  ahead.  The  presentation  on  the  echo  sounder  screen  is  a 
bird’s  eye  view  with  the  x-axis  pointing  starboard,  so  that  a  minus  sign 
must  be  applied  to  ux  in  formula  (13) 

It  should  be  noted,  that  the  amplitudes  of  the  three  signals  are  not  equal. 
The  three  transducer  sections  are  identical  in  geometrical  form  and 
tapering,  but  they  are  rotated  120  degrees  relative  to  each  other.  The 
geometrical  form  of  one  section  is  a  rhomb,  and  the  beam  pattern  is  not 
circular.  Therefore,  the  amplitude  of  the  three  signals  G1,  G2  and  G3 
should  not  be  used  individually. 

The  sum  of  the  three  echo  signals,  G1+G2+G3  represents  the  echo 
received  by  the  complete  transducer  area  and  is  used  for  the 
presentation  of  the  echo  amplitude  in  the  echogram,  and  in  the  target 
strength  calculation. 
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Beam  pattern  compensation 

The  main  lobe  of  the  transducer  ES  38-10  is  circular.  Thus  it  can  be 
expressed  as  a  function  of  vz  .  The  central  part  of  the  main  lobe, 
between  the  -6  dB  points,  can  as  a  good  approximation  be  expressed  by 
a  second  order  function 

B  ss  -0.126  Vz2  (14) 


where  vz  is  the  angle  from  the  z-axis  in  degrees 
B  is  one-way  beam  pattern  in  dB 


The  three  angles  from  s  to  the  co-ordinate  axes  must  obey 
cos2vx+  cos2vy+  cos2vz  -  1 
and  from  this 

sin2vz=  1-cos2vz  =  cos2vx+  cos2vy  =  sin2ux+sin2uy 

For  small  angles  with  si nu=u 

2  2  2 
Vz  =  U/+Uy 

B  =  -0.126  {uf+Uy2) 
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Transducer  design 

The  transducer  has  48  tonpilz  elements,  positioned  in  a  hexagonal 
pattern.  Each  tonpilz  element  is  circular  with  a  diameter  of  35  mm,  and 
the  centre  to  centre  distance  is  36.5  mm.  The  48  elements  are 
electrically  connected  as  three  sections  with  16  elements  in  each 
section.  The  transducer  centre,  origo,  is  a  point  between  the  three  inner 
elements.  The  circumference  of  the  48  elements  is  six-sided,  but  it  is  not 
regular.  The  sides  have  4  or  5  elements.  The  three  sections  are 
numbered  from  1  to  3: 

section  1,  starboard  aft 

section  2,  port  aft 

section  3,  forward 

The  transducer  elements  are  fixed  by  polyurethane  foam.  This  is  a  well- 
established  and  low  cost  technology. 


Transformers 

Three  transformers,  one  for  each  section,  are  placed  in  the  transducer, 
behind  the  tonpilz  elements.  The  transformers  tune  the  capacitive 
component  in  the  transducer  impedance,  and  transform  the  impedance  to 
75  ohms,  which  matches  the  impedance  of  the  transducer  cable.  In 
addition  the  transformer  is  used  for  amplitude  tapering  of  the  transducer 
elements  in  order  to  lower  the  side  lobes.  The  secondary  side  of  the 
transformer,  towards  the  elements,  has  four  taps  for  tapering  of  the 
elements.  The  voltage  levels  are  shown  on  the  drawing.  For  a  split-beam 
transducer'^  is  essential  to  avoid  that  a  large  fish  in  a  side  lobe  is  taken 
for  a  small  fish  in  the  main  lobe.  The  penalty  for  a  strong  tapering  is  a 
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wider  beam  and  limited  input  power  due  to  cavitation  at  the  transducer 
centre.  With  the  chosen  tapering  1500  W  may  be  applied  without 
cavitation  at  3  m  depth. 


Transducer  housing 


The  unit  with  the  tonpilz  elements  in  foam  is  completely  encapsulated  in 
a  polyurethane  mould,  which  makes  a  robust  and  watertight  housing.  The 
transducer  housing  is  circular  for  easy  installation  into  a  blister  on  the 
vessel  hull.  The  transducer  is  mounted  with  6  bolts  through  holes  in  the 
transducer  periphery.  The  holes  are  lined  with  sleeves  of  stainless  steel, 
so  that  the  bolts  can  be  tighten  securely,  metal  against  metal.  A  steel 
mounting  ring  has  to  be  welded  into  the  blister,  before  the  transducer  is 
installed. 
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Abstract 

Finite  element  (FE)  simulations  have  been  compared  with  classical  “Mason”  type  of  models  and  measurements 
for  electrical  response  functions  of  piezoelectric  disks.  It  is  shown  how  a  relative  converged  accuracy  of  the  FE 
simulations  better  than  the  relative  resolution  of  the  measurements  for  both  frequencies  and  response  values  for 
the  first  radial  modes  of  the  disk  can  be  obtained.  The  results  of  some  analyses  using  the  FE  model  with  adjusted 
constants  are  given,  particularly  for  the  lower  radial  modes  frequency  range.  Some  observations  on  the  effects  of 
frequency  dependent  dielectric  constants  are  also  given. 

1.  Introduction 

For  more  than  100  years  a  lot  of  work  has  been  done  on  trying  to  understand  and  simulate  the  behaviour  of 
piezoelectric  elements.  For  some  specific  resonance  modes  simplified  and  one-dimensional  analytical  models 
may  be  available,  such  as  the  much  used  “Mason”  type  of  models  for  thickness  extensional  (TE)  [1,2]  and  radial 
(R)  [2,3]  modes.  One  way  to  describe  the  various  material  loss  effects  is  through  the  use  of  complex  elastic, 
dielectric  and  piezoelectric  material  constants,  where  the  imaginary  parts  are  associated  with  the  losses  [4].  As 
these  solutions  are  valid  only  under  idealizing  assumptions,  the  models  can  not  describe  the  whole  mode 
spectrum  including  the  coupling  of  the  modes  for  practical  finite  dimensions  of  the  elements.  More  recently, 
numerical  finite  element  models  have  become  available  also  for  piezoelectric  structures  [5,6,7].  Such  models 
have  been  shown  to  provide  a  far  more  complete  and  accurate  description  of  piezoelectric  elements  and 
piezoelectric  transducer  constructions.  However,  the  use  of  such  models  for  accurate  simulations  is  still  limited 
due  to  limited  information  being  available  for  the  material  constants  involved,  and  also  due  to  limitations  in  the 
present  understanding  and  applications  of  how  the  different  constants  affects  the  element  vibrations. 

In  the  present  work  a  finite  element  (FE)  model  [7]  will  be  used  for  giving  some  examples  on  analyses  of 
electric  measurement  data  for  a  single  piezoelectric  disk  made  of  Morgan  Matroc  PZT-5A  material  [8].  The  disk 
is  denoted  “disk  21”,  and  the  diameter,  D,  is  40.10  mm  and  the  thickness,  T,  with  electrodes  is  1.989  mm,  see 
Fig.  1.  The  diameter  over  thickness  (D/T)  ratio  is  thus  20.1601.  Such  high  D/T  ratios  are  required  for  using  TE 
and  R  modes  in  methods  for  determining  material  constants  according  to  recommendations  given  in  e.g.  the 
IEEE  Std.  176  [2].  It  is  important  that  the 
numerical  uncertainty  of  the  FE  model 
calculations  will  be  well  below  the 
measurement  resolution  in  order  that  numerical 
artifacts  are  not  introduced  in  the  comparisons 
with  experimental  data.  This  aspect  is  given 
particular  attention  here,  and  some  results  are 
discussed  in  Secs.  3  and  4.  Further,  the  use  of 
the  FE  model  for  a  more  detailed  analysis  of 
the  electrical  measurements  of  disk  properties 
is  discussed  in  Sec.  5. 

An  example  of  typical  electrical  measurements  of  the  conductance  of  the  disk  21  is  shown  in  Fig.  2 
together  with  simulation  results  using  the  “Mason”  type  of  model  for  radial  modes.  Fig.  3  shows  the  same 
measurements  together  with  results  using  the  “Mason”  type  of  model  for  TE  modes.  In  Fig.  4  the  measurement 


PZT-5A  Electrodes 


Fig.  1.  The  disk  21  on  which  the  electrical 
measurements  and  simulations  have  been  performed. 


40 


data  are  compared  with  results  using  FE  simulations.  Note  that  logarithmic  scales  are  used  along  the  two  axes  in 
order  to  represent  the  large  span  in  frequency  and  conductance  values  involved.  For  the  first  radial  modes  a 
reasonable  agreement  with  the  radial  modes  theory  is  illustrated.  Also  a  reasonable  agreement  is  shown  with 
Mason  TE  model  at  the  first  TE  mode  in  Fig.  3.  The  limitations  of  the  two  models  used  in  Figs.  2  and  3  are 
clearly  seen  even  for  a  disk  with  such  a  high  D/T  ratio.  The  FE  model  is  shown  to  provide  at  least  a  more 
complete  qualitative  description  of  the  vibrations  of  the  disk.  The  quantitative  accuracy  in  the  simulations  may 
be  limited  through  the  limited  knowledge  available  concerning  the  material  constants,  and  the  limited  amount  of 
adjustments  being  done  here  for  the  results  shown.  This  is  one  of  the  main  motivations  for  this  type  of 
adjustments;  if  the  constants  are  not  given  with  the  desired  accuracy,  simulations  on  piezoelectric  elements  and 
transducer  constructions  will  be  inaccurate  and  often  useless  for  practical  purposes  [9],  Further  details  on  the 
models  and  the  data  used  will  be  found  in  Secs.  2  and  5. 


Fig.  2.  Comparison  between  measured  conductance 
and  calculations  using  Eq.  (1). 


Fig.  3.  Comparison  between  measured  conductance 
and  calculations  using  Eq.  (7). 


Fig.  4.  Comparison  between  measured  conductance  and  calculations  using  FE  simulations. 


2.  Theory  and  measurements 


The  “Mason”  type  of  model  used  to  describe  the  R  modes  in  Fig.  2  in  the  introduction  is  given  by  the  following 
expression  for  the  electrical  admittance,  Y  [10], 


Y(co)  =  G  +  iB  =  ico\ 


1- 


(l) 


In  Eq.  (1)  G  is  the  conductance  and  B  is  the  susceptance  of  the  element.  3(x)  =  xJ0(x)/JI(x)  is  called  the  “Ono’s” 
function,  a  =  D/2  is  the  radius  of  the  disk,  p  is  the  density  of  the  piezoelectric  material,  and  “  A  ”  is  used  for 
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denoting  complex  constants.  The  constants  used  in  this  model  are  derived  from  the 
FE  simulations,  see  Table  1,  by  [2], 
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Here  the  compliance  constant  matrix,  [s],  is  derived  from  the  stiffness  constant  matrix,  [c],  [s]  =  [c]'1  [2]. 

The  “Mason”  model  used  to  describe  the  TE  modes  in  Fig.  3,  is  given  by  the  following  expression  for  the 
electrical  impedance,  Z  [11], 


Z(co)  =  R  +  iX  = 


1  —  k, 


tan [a)T^p/cE  n) 


o)T^jp/cE  12 


i(om2i33 


(V 


In  Eq.  (7)  R  is  the  resistance  and  X  is  the  reactance  of  the  element.  The  constants  used  in  Eq.  (7)  are  derived 
from  the  set  of  constants  used  in  the  FE  simulations  as  given  in  Table  1,  by  [2], 
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All  the  constants  used  and  given  in  Table  1  have  the  same  definitions  as  in  [2],  except  that  complex  constants  are 
being  used  here  to  include  losses. 


The  program  FEMP  version  3.0  [7]  has  been  used  for  the  FE  simulations  in  the  present  work.  The 
program  can  calculate  eigenmodes,  eigenmode  resonance  spectra  [12],  element  displacements,  electrical 
response  functions,  and  also  sound  field  radiation  into  different  media  and  sensitivity  response  functions.  Front- 
and  backing  layers  in  more  complicated  transducer  structures,  can  be  included  [7].  Two  different  methods,  the 
direct  time  harmonic  method  and  the  mode  superposition  method,  for  calculating  the  response  functions  are 
implemented  in  the  FE  model  used  [7].  When  using  the  direct  method,  the  response  functions  are  calculated 
directly  from  the  FE  equations  using  matrix  manipulation.  For  the  mode  superposition  method  the  response 
functions  are  calculated  by  summing  the  contributions  to  the  response  for  all  eigenmodes  up  to  a  maximum 
frequency,  fniax  [7].  The  number  of  elements  per  shear  wavelength,  Ns,  of  sound  waves  in  the  radial  and  thickness 
directions  of  the  disk  at  given  frequencies  will  influence  on  the  accuracy  obtained  in  the  simulations.  For  the 
mode  superposition  method  f™*  also  affects  the  accuracy.  Various  number  of  elements  and  maximum  frequency 
has  been  chosen  in  the  different  examples  in  the  present  work.  Convergence  tests  of  the  convergence  accuracy 
using  the  FEMP  3.0  program  are  shown  in  Sec.  4. 

The  electric  measurements  on  the  disk  have  been  performed  using  a  Hewlet  Packard  4192A  impedance 
analyzer.  The  impedance  analyzer  has  been  controlled  from,  and  the  measurements  have  been  logged  on  a  PC 
via  a  standard  IEEE488-bus  and  a  Matlab  program.  A  vertical  holder  with  little  mechanical  loading  but  good 
electric  contact  has  been  used  to  hold  the  disk.  Different  frequency  resolutions  have  been  chosen  in  the  different 
examples  presented  here  (see  Sec.  5) 


3.  Accuracy  needed  for  the  FE  simulations 


In  Fig.  5  a  typical  measurement  series  close  to  the 
first  radial  mode  is  shown  to  give  an  example  on 
practical  measurement  resolution  both  in  frequency 
and  conductance  value  (similar  results  for  the 
susceptance  has  not  been  shown  here  for 
simplification).  A  measurement  resolution  of  1  Hz 
has  been  used  for  the  frequency.  The  series 
resonance  frequency  for  maximum  G  [2]  is  found  at 
50.531  kHz  with  a  relative  resolution  in  the 
determined  resonance  frequency  of  the  order  of  10 
ppm.  The  resolution  of  the  measurements  of  the 
conductance  value  is  0.01  mS  for  the  measurements 
shown,  which  represents  a  relative  resolution  of  the 
order  of  20  ppm.  (Note  that  in  the  present  work  the 


absolute  uncertainty  of  the  measurements  will  not  be 
considered  or  used). 
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Fig.  5.  High  resolution  measurements  of  the 
conductance  at  the  first  radial  serial  resonance. 


When  comparing  with  measurements  it  is  important  that  the  FE  calculations  can  be  done  at  a  higher 
accuracy  then  the  measurement  resolution,  in  order  to  ensure  that  significant  variations  in  the  responses  are  not 
due  to  calculation  errors.  Thus  the  considerations  of  the  measurement  data  show  that  an  accuracy  of  the  order  of 
1  ppm  is  of  interest  both  with  respect  to  resonance  frequencies  and  the  response  function  values.  How  such 
numerical  accuracies  can  be  achieved  in  the  FE  simulations  will  be  discussed  in  the  following  section. 


4.  Convergence 

Convergence  tests  for  the  FE  simulations  are  used  to  show  that  a  converged  accuracy  of  lppm  can  be  obtained 
for  the  first  radial  resonance  frequencies  in  both  frequency  and  value.  Such  convergence  tests  have  been 
performed  on  the  three  first  radial  resonance  frequencies  for  both  frequencies  and  conductance  values  and  also 
for  the  region  below  the  first  radial  serial  resonance.  The  converged  accuracy  is  found  for  both  the  direct  and  the 
mode  superposition  method.  Convergence  tests  are  performed  by  specifying  and  varying  Ns  for  both  methods.  In 
[7]  it  has  been  found  to  be  reasonable  to  choose  the  same  number,  Ns,  in  both  radial  and  thickness  directions.  For 
the  mode  superposition  method  convergence  tests  by  variations  of  fmax  has  been  performed  as  well. 

In  Fig.  6  fmax  for  the  mode  superposition  method  at  the  radial  modes  Rl,  R2  and  R3  has  been  chosen  to  be 
approximately  twice  the  respective  resonance  frequencies  fsRi  ®  50  kHz,  fsR2  ~  128  kHz  and  fsR3  «  201  kHz.  For 
the  direct  method  only  the  converged  accuracy  at  fSRi  has  been  shown.  Ns  at  fSRi,  fsR2,  or  fSR3  are  varied  from  Ns 
=  2  to  Ns  =  50  in  both  the  radial  and  the  thickness  directions  and  the  deviation  in  frequency  relative  to  the 
frequency  found  using  50  elements  per  shear  wavelength,  (fsRi.Ns  -  fsRi,Nso)/  f  sri,nso,  is  calculated.  1  ppm 
converged  accuracy  is  reached  at  Ns  approximately  12,  18  and  20  elements  pr  shear  wavelength,  respectively,  at 
the  three  resonances  for  the  mode  superposition  method.  For  the  direct  method,  1  ppm  converged  accuracy  is 
reached  at  Ns  =  9.  In'  these  examples  a  frequency  resolution  of  0.01  Hz  has  been  used  at  Rl  and  R2,  and  0.1  Hz 
at  R3.  This  resolution  makes  it  possible  to  detect  relative  differences  in  frequency  of  approximately  0.1  ppm. 

In  Fig.  7  the  deviation  in  conductance  values  at  the  resonance  frequencies  relative  to  the  conductance 
found  using  50  elements  per  shear  wavelength  at  the  resonance  frequencies  is  calculated  as  a  function  of  Ns. 
Convergence  tests  when  varying  fmax  for  the  mode  superposition  method  has  shown  that  the  conductance  value 
shows  no  sign  of  converging  even  for  =1.5  MHz.  This  means  that  for  the  example  in  Fig.  7  the  conductance 
value  is  converging  for  the  mode  superposition  method,  but  towards  the  wrong  value.  When  using  Ns  =  50  at  Rl, 
the  conductance  value  for  the  two  methods  deviate  with  152  ppm  when  fmax  =  100  kHz.  For  the  direct  method 
fmaX  is  of  no  influence,  and  the  conductance  is  converging  towards  the  correct  theoretical  conductance  value. 
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Fig.  6.  Deviation  in  frequency  for  the  three  first 
radial  resonances  relative  to  the  frequency found 
using  50  elements  per  shear  wavelength. 


Fig.  7.  Deviation  in  conductance  values  for  the  three 
first  radial  resonances  relative  to  the  values  found 
using  50  elements  per  shear  wavelength. 


When  using  the  mode  superposition  method  this  is  only  a  problem  for  the  conductance  values  and  not  for  the 
frequencies.  The  convergence  is  seen  to  be  a  little  slower  for  the  conductance  in  Fig.  7  than  for  the  resonance 
frequencies  in  Fig.  6,  and  a  1  ppm  converged  accuracy  is  reached  at  Ns  of  approximately  16,  22  and  28  elements 
per  shear  wavelength,  respectively,  at  the  three  resonances  for  the  mode  superposition  method.  For  the  direct 
method,  1  ppm  converged  accuracy  is  reached  at  Ns  «■  13. 

With  the  exception  of  the  convergence  for  the  conductance  as  a  function  of  f^  for  the  mode 
superposition  method  all  tests  show  1  ppm  converged  accuracy  with  less  than  28  elements  per  wavelength  at  the 
frequencies  tested.  It  is  of  no  problem  to  calculate  with  this  number  of  elements  at  these  frequencies  using  a  PC 
with  a  550  MHz  processor  and  256  MB  of  internal  memory.  1  ppm  converged  accuracy  is  therefore  fully 
achievable  both  for  frequency  and  conductance  values  at  the  first  radial  modes  and  for  the  conductance  values  at 
the  lower  frequencies.  Another  question  is  of  course  whether  the  FE  model  simulates  the  response  functions 
physically  correctly  to  such  an  accuracy.  The  FE  model  used  in  this  work  has  been  compared  with  other 
commercial  FE  models  with  a  good  agreement  [13].  But  only  detailed  comparisons  with  measurements  at  a 
sufficiently  high  accuracy  can  tell  how  accurate  the  FE  model  really  is.  Such  an  endeavor  has  not  been  attempted 
here. 


5.  Analyses 

In  the  FE  simulations  shown  in  Fig.  4  the  direct  method  with  Ns  =  4  at  2  MHz  was  used  to  obtain  the 
conductance  response.  To  get  the  agreement  shown  it  was  necessary  to  adjust  die  constants  used  in  the  FE 
simulations.  It  should  be  noted  that  all  the  adjustments  done  in  this  section  are  partly  of  a  “serendipity  type”,  as  a 
full  understanding  of  how  the  different  constants  affect  the  response  functions  at  a  high  accuracy  is  lacking.  Thus 
the  values  of  the  constants  arrived  at  through  die  adjustments  are  not  necessarily  the  physically  correct  values. 
Variations  were  done  using  as  a  starting  point  a  FE  set  of  real  constants  found  in  Chapter  5  in  [14]  to  approach  a 
radial  mode  set  of  complex  constants  derived  by  Eqs.  (2)-(6),  as  also  found  in  [14],  Table  6.5.  Thus  the  approach 
is  to  first  obtain  an  adjustment  in  the  low  frequency  radial  modes  range  based  on  adjusting  the  four  planar 

constants  £33 ,  c/J ,  S~p  and  kp  in  the  radial  modes  model  in  Eq.  (1).  In  this  process  cfj  was  adjusted  to  get 
matching  of  cp ,  £3i  was  adjusted  to  get  matching  of  ££ ,  e3,  was  adjusted  to  get  matching  of  kp ,  and  cf2 
was  adjusted  to  get  matching  of  ar .  The  piezoelectric  constants,  e,  were  all  assumed  to  be  real  before  any 
adjustments  were  performed,  but  a  complex  e3l  has  been  chosen  here  to  adjust  kp ,  while  ey  and  e/s  were 
kept  real.  To  get  an  even  better  agreement  with  the  measurements  for  the  first  R-modes  as  seen  in  Figs.  8  and  9, 
more  accurate  adjustments  of  the  real  part  of  CtI  and  the  imaginary  part  of  c^2  were  performed.  It  is  these 
adjusted  constants  and  the  constants  derived  from  these  which  are  given  in  Table  1  and  which  have  been  used  in 
the  simulations  for  all  of  the  Figs.  2,  3,  4,  8  and  9.  It  must  be  stressed  that  these  adjustments  only  show  one 
possible  way  to  get  better  agreement  between  the  simulations  and  the  measurements,  and  is  not  necessarily  the 
physically  correct  choice,  as  noted  above. 

Table  1.  The  material  constants  needed  in  the  FE-simulation,  R-mode  simulations  and  TE-mode  simulations.  The 
PZT-5A  column  is  constant  values  given  by  Morgan  Matroc  [8J.  The  last  column  shows  the  adjusted  constant 
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values  that  are  used  in  the  simulations.  The  R  and  TE  -mode  constants  are  derived  from  the  FE-set  of  constants 
by  Eqs.  (2)-(6)  and  Eqs.  (8)  and  (9). 


Symbol 

Unit 

PZT-5A 

Adjusted  constants 

Symbol 

PZT-5A 

Ad  j.  constants 

FE-const. 

R-modes 

Derived  values 

P 

Hgg 

7750 

7740 

- 

- 

1040x(l-i/93) 

^11  /£o 

916 

916x(l-i/50) 

| 

-0.60 

-0.5834x(l-i/122) 

- 

830 

825x(l-i/92) 

IM 

gmn 

- 

7.25308x(l+i/98) 

!■ 

-5.4 

-5.59x(l-i/44) 

dp 

- 

0.3541x(l-i/485) 

e33 

II 

15.8 

14.79 

| 

Derived  values 

ei5 

It 

12.3 

12.3 

£*3/£0 

- 

830 

825x(l-i/92) 

gum 

12.1 

10.9835x(l+i/98) 

K 

- 

0.48 

0.454x(l+i/3788) 

4 

it 

7.54 

7.30x(l+i/106.732) 

ci 

14.7 

14.506x(l+i/97) 

CB 

‘•'13 

II 

7.52 

7.38x(l+i/98) 

cE 

‘■'33 

it 

11.1 

11.51x(l+i/98) 

cE 

U44 

It 

2.11 

2.11x(l+i/98) 

tan  d 

- 

0.02 

- 

Qm 

- 

75 

- 

More  detailed  comparisons  around  R1  in  Figs.  8  and  9  show  that  using  the  derived  sets  of  constants,  the 
radial  model  results  are  shifted  a  little  above  the  FE  simulations  in  frequency  as  is  also  to  be  expected.  In  the 
Mason  model  an  infinitely  large  and  thin  plate  is  assumed,  and  deviations  between  the  FE  simulations  and  the 
results  from  the  radial  modes  model  can  be  seen  even  for  D/T  =  20.  This  deviation  will  be  larger  for  the  higher 
R-modes  as  seen  in  Fig.  3.  A  Rayleigh  type  of  correction  [15]  could  be  used  to  reduce  this  problem,  but  the  IEEE 
176-1987  standard  [2]  does  not  recommend  the  use  of  such  a  Rayleigh  correction  when  the  main  intention  of  the 
measurements  is  the  determination  of  material  constants.  For  the  FE  simulations  in  Figs.  8  and  9  the  direct 
method  for  calculating  the  conductance  with  Ns  =  20  at  50  kFIz  has  been  used.  For  the  measured  data  in  Figs.  8 
and  9  a  limited  frequency  resolution  of  20  Hz  has  been  used.  A  much  better  frequency  resolution,  such  as  the  1 
Hz  resolution  used  for  Fig.  5  could  have  been  used  to  obtain  a  smoother  and  more  accurate  conductance  curve. 
(Less  accurate  results  using  a  resolution  of  100  Hz  are  given  in  Figs.  2,  3,  and  4).  It  should  be  mentioned  that  if 
these  comparisons  had  been  done  using  the  constant  values  provided  by  Morgan  Matroc  [8],  much  larger 
deviations  from  the  measurements  would  have  been  observed.  It  should  also  be  noted  that  no  possible  effects 
from  finite  electrode  thickness  have  been  included  in  the  present  analyses. 


Fig.  8.  Comparisons  between  measured  conductance  (dotted  line),  simulations  using  the  Mason  radial  mode 
model,  and  simulations  using  FE. 


Fig.  9 .  Comparisons  between  measured  conductance  Fig.  10 .  Measured  capacitance  and  simulations  with 
(dotted  line),  simulations  using  the  Mason  radial  mode  and  without  frequency  dependence  using  the  Mason 

model,  and  simulations  using  FE.  radial  mode  model. 


If  the  response  is  plotted  on  the  form  of  capacitance  (B/ro),  the  effects  of  a  frequency  dependent 
dielectric  material  constant  can  be  seen.  For  the  Mason  model  it  is  easy  to  implement  such  a  frequency 
dependence.  In  Fig.  10  the  measured  and  calculated  capacitance  is  plotted  without  and  with  a  frequency 
dependent  et33.  et33  has  been  fitted  at  500  Hz  and  reduced  by  2.05%  per  frequency  decade  to  fit  the  trend  in  the 
experimental  data.  Morgan  Matroc  [8]  gives  a  frequency  dependence  of  eT33  to  be  -2,4%  per  frequency  decade. 
This  effect  of  frequency  dependence  is  shown  for  the  low  frequency  range  in  Fig.  10,  but  is  also  expected  to  give 
further  deviations  between  measurements  and  simulations  for  higher  frequencies  as  well.  Until  now,  no  FE 
simulations  including  such  frequency  dependent  constants  has  been  performed. 

6.  Conclusions 

The  FE  simulations  of  the  electrical  response  functions  for  piezoelectric  disks  are  seen  to  give  a  much  better 
overall  agreement  with  the  measurements  compared  to  the  classical  “Mason”  type  of  models.  The  FE  model  has 
been  shown  to  be  capable  of  calculations  at  a  higher  converged  accuracy  than  the  resolution  of  the  measurements 
for  the  frequency  range  of  the  first  radial  modes  and  below.  This  enables  very  accurate  adjustments  of  the 
material  constants  that  describes  the  piezoelectric  disk  by  comparison  between  measurements  and  simulations  of 
the  electrical  response  functions.  More  extensive  FE  analyses  of  the  behaviour  of  piezoelectric  elements,  will 
benefit  from  further  detailed  work  on  studying  the  effects  of  the  different  specific  material  constants  on  the 
element  vibrations. 
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Ultrasonic  flowmeters  are  increasingly  being  used  in  the  fiscal  metering  of  dry  natural 
gas.  Typically  such  meters  are  in  the  size  range  of  6  to  40”  diameter  and  can  offer 
accuracy  of  0.5  %  or  better.  Currently  there  is  much  interest  in  the  development  of 
this  technology  for  use  in  the  metering  of  wet  natural  gas.  This  presents  several 
challenges  if  good  accuracy  is  to  be  maintained.  This  abstract  describes  one  aspect  of 
the  development  of  ultrasonic  wet  gas  metering,  that  is  the  measurement  of  ultrasonic 
scattering  from  liquid  droplets  and  mist  carried  in  the  gas.  This  work  is  part  of  a  3- 
year  project  jointly  funded  by  Kongsberg  Offshore,  Statoil,  Norsk  Hydro,  Phillips 
Norge  and  the  Norwegian  Research  Council  (NFR). 

A  considerable  amount  of  literature  exists  on  the  many  theoretical  models  for 
ultrasonic  scattering  from  single  scatterers  and  clouds  of  scatterers.  These  models 
have  been  tested  widely  against  experimental  measurements  in  the  fields  of 
underwater  acoustics  and  medical  ultrasound  where  the  continuous  medium  is  liquid 
and  the  dispersed  medium  is  a  gas  or  a  solid.  There  have  also  been  comparisons  of 
such  theoretical  models  in  emulsions.  Although  there  is  no  reason  to  suspect  that  the 
theoretical  models  are  inapplicable  to  the  gas-continuous  case,  there  is  a  dearth  of 
publications  in  the  open  literature  on  measurements  of  ultrasonic  scattering  from 
liquid  droplets  or  mist  in  suspension  in  a  gas. 

In  this  study  we  have  attempted  measurements  of  ultrasonic  scattering  (side  scatter 
and  backscatter)  on  a  wide  range  of  droplet  sizes  (1  micron  to  8  mm  diameter)  and 
liquid  volume  fractions  (10  -  1155  parts  per  million)  for  ultrasonic  frequencies  in  the 
range  50  to  300  kHz. 

Measurements  of  side  scatter  and  back  scatter  made  on  single  plastic  spherical  beads 
(diameters  in  the  range  3  to  8  mm)  showed  good  agreement  with  theoretical 
predictions  (within  the  measurement  uncertainty  of  about  2  dB).  Measurements  of 
back  scatter  from  a  liquid  mist  (droplet  diameters  33  -  52  micron)  gave  scattering 
levels  some  15  to  20  higher  than  predicted.  It  is  speculated  that  the  discrepancy  could 
be  due  to  interface  effects  at  the  boundaries  of  the  mist  filled  region  and/or  scattering 
from  turbulence  caused  by  the  spray  nozzles  used  to  generate  the  liquid  mist.  One  of 
the  main  difficulties  of  such  measurements  is  generation  of  a  well  characterised  mist, 
the  two  main  parameters  of  interest  being  the  droplet  diameters  and  the  volume 
fraction  occupied.  In  an  effort  to  create  a  well  defined  scattering  medium  an  array  of 
spherical  plastic  beads  (4  mm  diam)  suspended  on  fine  wires  was  built.  The  beads 
were  placed  in  a  random  pattern  and  the  suspension  wires  were  small  enough  that 
they  didn't  make  a  significant  contribution  to  the  scattered  field.  Back  scattering 
measurements  made  on  this  random  array  showed  better  agreement  with  theory 
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(within  about  5  dB)  but  the  difference  was  still  outside  of  the  calculated  range  of  the 
expanded  total  uncertainty. 

It  is  important  to  obtain  good  experimental  measurements  in  order  to  verify  the 
candidate  theoretical  models.  Such  theoretical  models  will  then  have  the  advantage 
that  they  offer  greater  flexibility  and  are  easy  to  apply  under  conditions  where  good 
experimental  measurements  would  be  difficult  to  obtain,  e.g.  at  high  pressures  and 
temperatures. 


A.  C.  Baker,  15  March  2001. 
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Abstract 

A  one-dimensional  description  of  reciprocal  ultrasound  systems  is  described. 
The  model  allows  near-  and  far-fleld  characteristics  to  be  evaluated  under  various 
excitation  and  electrical  loading  conditions.  Numerical  results  are  obtained  for 
reciprocal  transducer  systems  of  interest  for  flow  meter  district-heating  applications 
using  transducers  with  half-wavelength  thick  stainless  steel  mstchtwg  layers.  Good 
agreement  with  measurements  is  found. 


1  The  piezoelectric  plate  transducer  operated  in  thickness 
mode 

The  following  discussion  is  restricted  to  piezoelectric  transducers  in  the  form  of  a  thin  cylinder.  We 
shall  assume  that  the  piezoelectric  transducer  is  operated  in  thickness  mode.  Under  these  simplifying 
circumstances,  the  general  3D  piezoelectric  equations  degenerate  into: 


D  =  dT  +  eT  E 

(1) 

E  =  pTD  -  gT 

(2) 

D  =  esE  +  eS 

(3) 

E  =  fdsD  -  hS 

(4) 

S  =  dE  +  seT 

(5) 

T  =  -eE  +  cES 

(6) 

S  =  gD  +  sdT 

(7) 

T  =  -hD  +  cDS. 

(8) 

It  should  be  mentioned  that  two  independent  relations,  e.g.,  Eqs.  (3)  and  (8),  suffice  to  describe  the 
one-dimensional  problem.  The  rest,  Eqs.  (l)-(2)  and  (4)-(7),  are  just  alternative  representations  of 
Eqs.  (3)  and  (8)  in  different  variables. 


2  Dynamic  operation 


Newton’s  Second  Law,  Poisson’s  Law,  and  the  definition  of  strain  completely  determine  the  temporal 
and  spatial  behavior  of  the  resulting  strains,  stresses,  velocities,  currents,  and  voltages  in  the  trans¬ 
mitter  and  receiver  transducers  when  invoking  the  appropriate  mechanical  boundary  conditions  and 
subsidiary  electromechanical  conditions  (e.g.,  the  coupling  between  applied  voltage  and  electric  dis¬ 
placement,  transducer  electrode  velocities).  In  this  section,  the  one- dimensional  modeling  equations 
accounting  for  the  dynamic  response  of  a  multilayer  transmitter-receiver  piezoelectric  transducer 
setup  will  be  derived. 

Newton’s  Second  Law  applied  to  an  infinitesimal  slab  of  width  dz  can  be  stated  as  [4]: 

pAdz^  =  (Fz  -  Fz+dz)  =  -  AT,  +  ATz+dz  =  A^dz,  (9) 

where  u  is  the  particle  velocity,  p  is  the  mass  density,  A  is  the  cross-sectional  area,  and  F,  is  the 
force  on  the  slab  at  position  z.  Eq.  (9)  is  equivalent  to: 


du_dT 
P  dt  dz’ 


(10) 


and  from  the  definition  of  strain  in  the  one-dimensional  case:  S  —  it  is  found  that 


du  9  (ft )  d  d(  _  dS 
dz~  dz  ~  dt  dz~  dt' 


where  <  is  the  z-coordinate  displacement.  Equations  (10-11)  are  the  basic  (one-dimensional)  equa¬ 
tions  describing  the  dynamics  of  a  piezoelectric  material  and,  of  course,  an  ordinary  (non-piezoelectric) 
material.  Accordingly,  these  equations  apply  to  a  multilayer  transducer  consisting  of  a  piezocerainic 
material  as  well  as  ordinary  materials.  It  is  convenient  to  cast  Eqs.(lO-ll)  in  a  different  form  by  use 
of  the  Poisson  Equation: 


and  Eq.(8)  implying: 


9D  n 

~di  ~  Pfree  ~  °’ 

(12) 

du  DdS 

P  dt  °  dz 

(13) 

du  dS 

(14) 

dz  dt 
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In  Eq.(12),  pjree  denotes  the  external  (or  free)  charge  density,  which  is  zero  everywhere  within  the 
multilayer  transducer  except  on  the  electrodes.  Eqs.  (13)  and  (14)  therefore  apply  to  the  piezoccramic 
layer  as  well  acoustic  coupling  layers  and  matching  layers. 

The  voltage  of  the  transmitting  voltage  source  is  related  to  the  electric  field  between  the  electrode 
plates  as  follows: 


V(t)  =  ZelI+f  E{z,t)dz ,  (15) 

J  ZL 

where  zi  and  zj  are  the  left  and  right  ceramic  plate  positions  of  the  transmitter  multilayer  structure, 
respectively  (see  Fig.l),  Zel  represents  any  complex  electrical  impedance  in  series  with  the  pulse 
generator  and  the  transmitting  transducer,  and  /  is  the  current.  The  impedance  Ze i  may  account 
for  internal  electrical  impedances  in  the  pulse  generator  as  well  as  external  ohmic  resistors  and/or 
inductive  or  capacitive  impedances.  The  loads  on  the  left  and  right  boundaries  of  the  transmitter 
are  considered  to  be  air  and  water,  respectively. 

Differentiating  Eq.(15)  with  respect  to  time  and  making  use  of  Eqs.  (4),  (12),  and  (14)  leads  to 
(keeping  in  mind  the  relation  between  current  and  electric  displacement:  I  —  A^-): 

8V  d2D  hdD  Lr  .  v  „ 

~ai  =  AZel  ~W  +  7s-Q[-h[<^)-^L)\,  (16) 

where  A  denotes  the  cross-sectional  area  of  the  piezoccramic  material.  Eq.  (16)  is  the  electrome¬ 
chanical  subsidiary  condition  needed  for  the  transmitter.  An  equation  similar  to  Eq.  (16)  can  be 
written  down  for  the  receiver;  refer  to  Eq.  (34)  and  Fig.2. 

The  set  of  equations:  (13),  (14),  (16)  and  continuity  of  the  (normal)  velocity  and  pressure  between 
layers  and  surrounding  media  completely  specify  the  problem  to  be  solved.  However,  before  stating 
these  conditions  explicitly,  some  further  discussions  will  be  carried  out.  Firstly,  notice  that  Eqs.(13) 
and  (14)  can  be  augmented  in  S  and  u,  respectively,  to  yield: 


dlu  .,  d2u 

dF~Vadz*~° 


&2s 
at 2 


_  _  _  2  d'2S  _ 

~Va~azi 


(17) 

(18) 


where  va  =  y  ^  is  the  layer-dependent  sound  velocity,  i.e.,  it  depends  on  the  position  coordinate  z 
for  a  multilayer  structure. 
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Figure  1:  Schematic  transducer  geometry,  transmitter. 


In  the  present  work,  the  multilayer  structure  is  assumed  to  consist  of  a  piezoccramic  material 
(thickness  lpz),  an  acoustic  coupling  layer  (e.g.,  grease,  thickness  lac),  and  a  matching  layer  (thickness 
Ima)-  To  simplify  matters,  we  have  disregarded  the  thicknesses  of  the  plate  electrode  layers  in  this 
example  but  the  arguments  involved  can  easily  be  extended  to  any  multilayer  transducer  consisting 
of  an  arbitrary  number  of  layers  and  thicknesses.  In  fact,  it  has  been  verified  numerically  [5]  that  a 
few  microns  thick  electrode  layer  has  very  little  or  no  effect  on  results. 

2.1  Solving  transducer  dynamic  equations  by  use  of  Fourier  analysis 

The  formulation  of  boundary  conditions  are  conveniently  done  for  each  frequency  component  at  a 
time.  Applying  Fourier  analysis  to  the  generator  input  voltage  is  the  standard  method  which  will  be 
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used  here.  A  general  excitation  voltage  V(t),  provided  by  the  pulse  generator,  can  be  decomposed 
into  its  Fourier  components  V(u)  according  to: 


1  /■“ 

V(t)  =  -y —  /  V {u)  exp(-iuit)duj 

V27T  J-oo 

(19) 

1  f°° 

V(u)  =  -1=  V(t)oxp(iut)dt. 

V2x  J-oo 

(20) 

Solving  for  each  frequency  component  (■^^^(w)  exp (—iut)d(x>)  of  V  (t)  at 
ential  equations  (17)-(18)  reduce  to  ordinary  differential  equations: 

a  time,  the  partial  differ- 

(21) 

(22) 

with  the  shnple  solutions: 

u  =  uA  exp  (—ikz  —  iut)  +  ub  exp(ikz  —  iut) 

(23) 

S  -  Sa  exp  (-ikz  -  iut)  +  Sb  exp  (ikz  -  iut). 

(24) 

Here,  the  wave  number  k  =  —  has  been  introduced,  and  ua,  ub  ( Sa ,  Sb)  are  specified  by  the 
boundary  conditions.  It  follows  from  Eqs.(23)-(24)  that  the  plane  waves  associated  with  uA,  SA 
propagate  backwards  (hi  the  —z  direction)  whereas  ub,  Sb  propagate  forwards  (in  the  +z  direction). 
It  is  convenient  to  write  down  similar  expressions  for  the  stress  and  velocity  functions  in  terms  of 
one  set  of  coefficients  only,  say  Sa,  Sb-  According  to  Eqs.  (8),  (13),  and  (23)-(24),  the  sought 
expressions  become 

T  =  cP  [SA  exp (—ikz  —  iut)  +  Sb  exp (ikz  —  iut)]  —  hD  (25) 

u  =  -~A  exp  {-ikz  -  iut)  -  °  exp  {ikz  -  iut),  (26) 

Za  "O 

where  Za  -  pva  (note  that  for  plane  waves:  Za  =  pva  =  |£|).  Furthermore,  under  monofrcquency 
conditions,  Eq.(16)  simplifies  to: 

—iuiV  =  -AZelJ2D  -  -  h  [u(zi)  -  u(zL)\ .  (27) 


Parameter 

Value 

Unit 

u0/(2tt)  (driving  frequency) 

1-10B 

Hz 

N  (number  of  periods) 

8 

V0  (applied  voltage  amplitude) 

10.0 

V 

Thickness  (piczocerainic  material) 

1.96 

mm 

Transducer  area 

380 

min2 

Thickness  (stainless  steel) 

2.81 

min 

h  (piezocerainic  material) 

1.42  10® 

F2V/m3 

es  (piezoceramic  material) 

1440eo 

F/m 

cD  (piezocerainic  material) 

1.19 -1011 

Pa 

cD  (stainless  steel) 

2.53  •  1011 

Pa 

cP  (water) 

2.25  •  109 

Pa 

cD  (air) 

1.24  •  105 

Pa 

p  (piezoceramic  material) 

7750 

kg/m3 

p  (stainless  steel) 

8000 

kg/m3 

p  (water) 

1000 

kg/m3 

p  (air) 

1.29 

kg/m3 

Table  1:  Table  1:  Material  parameters  and  characteristic  dimensions. 
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Employing  the  two  equations  (25)- (26)  and  imposing  continuity  of  the  (normal)  velocity  and 
pressure  between  layers/surrounding  media  supplemented  by  the  condition  given  in  Eq.(27)  yield  nine 
equations  in  nine  unknowns  for  the  multilayer  transmitter  structure  depicted  in  Fig.l.  The  unknowns 
arc:  D,TL,TR,  SAi,SBni  =  1,2,3  where  TL  and  TR  refer  to  the  stress  amplitude  associated  with 
the  backward-propagating  wave  L  (in  air)  and  the  forward-propagating  stress  wave  R  (in  water), 
respectively.  The  explicit  form  of  the  nine  equations  read: 


Ti  —  cf  (.SU,  +  Sbi)  —  hi  D 

—L  -  _  £gi  \ 

ZaL  1  \  Za i  Za i  / 

cf +  Sb2)  =  cf  (SAl  exp(-ife1;1)  +  5b,  exp(ifc1;1)) 

-  hiD 

cd  ( 5^  _  SbA  =cp  ( sAi  cxp(-ikih)  _  SBl  exp{ikih)\ 

\Za2  Za2 )  1  \  Za I  Za  1  / 

cf  isA,  +  Sb3)  -  cf  (SA,  cxp(-ik2l2)  +  Sb2  cxp(ik2l2)) 

,p  ( SA3  Sb3  \  exp (~ik2l2)  Sb2  exp(ik2l2)\ 

^  U-3  ZaJ-^{  J 

Tn  =  cf  (S„3  exp(-ik3l3)  +  Sb3  exp {ik3l3)) 

=  P  ( Sa3  cxp(~2^^)  _  Sb3  exp (ik3l3)\ 

ZaR  \  Za  3  Za  3  / 

-,W  =  -AZ,1U‘D  -  - H,  [cf  (|j  -  fe)  ~  *  (fe  -  If)]  ’  (28) 

where  cf ,  Zam,  km,  hmi  cf  are  the  bulk  modulus  under  constant  electric  displacement,  acoustic 
impedance,  wave  vector,  piezoelectric  h  constant,  permittivity  under  constant  strain  of  layer  m, 
respectively,  and  ZaR  (ZaL)  is  the  acoustic  impedance  of  the  medium  located  right  (left)  to  the 
transducer.  Keep  in  mind  that  layer  1,  layer  2,  and  layer  3  in  the  case  of  the  transmitter  arc  the 
piezoceramic  layer,  the  acoustic  coupling  layer,  and  the  matching  layer,  respectively,  so  that,  e.g., 
li,  l2,  and  l3  are  equal  to  lp:,  lac,  and  lma,  respectively. 

Solving  Eqs.(28)  for  TR,  the  right  aperture  velocity  amplitude  uR  can  be  expressed  as: 


(29) 


and  insertion  into  the  Rayleigh  equation  gives  the  incoming  pressure  wave  at  the  receiving  trans¬ 
ducer  [6], [7]: 


Prec(x',W,t)  : 


UR(t  -  $  -  Z\ /van) 


assuming  ultrasound  transmission  between  the  two  transducers  occurs  in  a  water  bath  and  under 
zero-flow  conditions.  In  Eq.(30),  Atra,  vaR,  x,  and  x'  arc  the  transmitter  aperture  area,  sound  speed 
in  the  medium  to  the  right  (water),  the  position  coordinate  on  the  transmitter  aperture  (variable 
during  integration),  and  the  position  coordinate  on  the  receiver  aperture.  Notice,  that  the  axial 
distance  between  the  two  transducers,  i.e.,  the  z-coordinatc  of  x'-x,  is  a  constant  during  integration 
due  to  the  assumption  that  the  two  transducers  point  facc-to-facc.  Since  this  is  a  one-dimensional 
model,  the  incoming  pressure  on  the  receiver  aperture  (being  an  input  parameter  to  the  receiver)  is 
calculated  as  an  average  pressure  over  the  receiver  aperture: 


=  J-f 

Arec  J  Arte 


Trecix  ,  CU,  t)(Z.Arec, 


where  Arec  denotes  the  receiver  aperture  area.  Following  this  procedure  for  determining  incoming 
receiver  pressures  allows  receiver  voltage  signals  to  be  calculated  whether  the  receiver  is  located  in 
the  near-field  or  in  the  far-ficld  zones  of  the  transmitter.  In  a  previous  paper,  however,  the  receiver 
was  assumed  to  be  located  in  the  far-field  zone  of  the  transmitter  [5]. 
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In  the  case  where  transmission  takes  place  in,  e.g.,  a  flow  meter  in  a  moving  medium,  the  analysis 
leading  to  the  received  pressure  wave  becomes  much  more  involved.  When  comparing  stresses  instead 
of  pressures,  the  input  signal  to  the  receiver  is: 

Ti  =  Trec(w,  t)  =  -Precis  t).  (32) 


Ima  lac  Ipz 


ZL=Z0  Z1  22  ZR=Z3 


Layer  1:  Layer  2:  Layer  3: 

Matching  Acoustic  Piezoceramic 

layer  coupling  material 

layer 


Figure  2:  Schematic  transducer  geometry,  receiver. 


The  analogue  set  of  boundary  expressions  for  the  receiver,  with  Fig.2  as  a  reference,  can  now  be 
formulated  as: 

Ti  +  Tb  =  c?(SCl+SDl) 

+  p(Scl_SdA 

Zal  ZaL  1  \  Zal  Zal) 
cf  (Sc2  +  So2)  =  cf  (SCl  exp(-ifcili)  +  SDl  exp(iMi)) 
n  ( Sc,  Sd 2  \  n  ( Sch  exp(-ifciZi)  SDl  exp(?7,:i/1)  '\ 

C2  Ua2  Za2)-Cl{  Zal  Zal  ) 

cf(Sc3  +  SD3)  -  h3D  =  cf  (Sc,  exp(-ik-2h)  +  Sd,  exp (ik2h)) 
n  ( Sc,  So,  \  D  (Sc,  exp(-ik2h)  Sn2  cxp(ik2h)\ 

Cs  U-3  ZaJ-C2{  Za2  Za:>  ) 


Tf  =  cf  (Sc3  exp(—ik3l3)  +  Sd3  exp(ik3l3))  —  h3D 
Tf  _  p  ( Sc3exp{-ik3l3 )  _  5g3  exp(ik3l3)  \ 

ZaR  3  \  Za  3  -^a3  / 


where  Tb  and  Tf  are  the  reflected  stress  amplitude  at  the  left  aperture  and  the  forward  propagating 
stress  amplitude  at  the  right  aperture,  respectively,  and  Ze2  represents  the  scries  electrical  impedance 
in  the  receiver  circuit.  Note  that  in  the  case  of  the  receiver,  layer  1,  layer  2,  and  layer  3  arc  the 
matching  layer,  the  acoustic  coupling  layer,  and  the  piczoceramic  layer,  respectively  so  that  Zi ,  l2, 
and  l3  are  equal  to  lma,  lac,  and  lpz,  respectively,  etc.  The  time-dependent  analogue  to  the  last  of 
Eqs.(33)  is: 

0  -  7lZe2  |  ™  -  h  [«(**)  -  «(*»)] .  (34) 

Solving  Eqs.(33)  for  the  electric  displacement  D,  the  induced  voltage  Vf.ec(u,t)du  may  be  found: 

Vr'ec(w,  t)du  =  Z2I  =  -iZ2AuD,  (35) 

and  summing  up  contributions  from  all  frequency  components  associated  with  the  applied  voltage 
pulse  gives  the  real-time  induced  voltage  Vrcc{t): 

Vrec(t)  =  [  Vtec{w,t)du  -  —4=  [  Vrec(u,t)oxp(-iut)dw,  (36) 

J — oo  yj 2iK  J — oo 

where  Vrec(uj,t)  in  the  latter  equality  is  the  Fourier  transform  analogue  to  Vrec(t)  (compare  F(w) 
with  V(t)  in  Eqs.  (19)-(20)).  This  concludes  the  derivation  of  the  received  voltage  response  from  a 
known  applied  voltage  signal. 
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Figure  3:  Pulse  generator  voltage  (transmitter  circuit);  b)  Fourier  transform  of  pulse  generator 
voltage,  (N  =  8,  V0  -  10  V,  u0/{2ir)  =  1  MHz). 


3  Transmitter  excitation 


Consider  the  case  where  a  sinusoidal  8-pcriod  voltage  pulse  V'(<)  with  amplitude  F0  =  10  V  and 
frequency  w0/(2?r)  =  1  MHz  is  applied  to  the  piezoccrarnic  electrodes  equivalent  to: 

V(t)  =  V0  sin  w0t,  if  \t\  <  — 

ou'o 

V(t)  =  0,  if  |t|>-.  (37) 

Wo 

Deconvoluting  V (i)  into  its  Fourier  components  V(w)  according  to  Eq.  (20)  yields: 


F(w) 


iJlv° 


sin[(w0  -  w)(87r/w0)] 


2(w0  -  w) 

sin[(w0  +  w)(87t/wq)]] 


2(wq  +  w) 


(38) 


The  applied  real-time  voltage  is  shown  in  Fig.  3a  and  the  corresponding  Fourier  transform  is 
depicted  in  Fig.  3b.  The  resonance  is,  of  course,  located  at  the  driving  frequency:  1  MHz  (or  -1 
MHz)  and  falls  off  in  an  oscillating  manner  away  from  the  resonance. 

The  design  of  the  reciprocal  transducer  is  schematically  illustrated  in  Figs.  1-2,  and  the  various 
material  parameters  and  layer  thicknesses  are  given  in  Table  1.  Note  that  the  thicknesses  of  the 
piczoceramic  material  and  stainless  steel  are  chosen  to  be  exactly  a  half  wavelength  at  1  MHz  in  the 
present  work,  i.e.,  t£  =  /pi/(u0/2n);i  =  1,3.  These  values  are  slightly  different  from  those 

used  in  [8]. 


Figure  4:  Voltage  across  the  transmitter  electrodes  for  the  multilayer  structure  shown  in  Fig.  (1) 
where  the  acoustic  coupling  layer  thickness  is  zero,  a)  Calculated  voltage  response  b)  Measured 
voltage  response  under  similar  conditions.  Zi  =  75  Ohm,  Z2  =  75  Ohm. 


Measurements  were  carried  out  using  a  LcCroy  9400  DUAL  digital  oscilloscope  to  aquire  voltage 
signals  across  the  receiver  electrodes.  The  transmitter  was  electrically  excited  using  a  Philips  PM 
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138  function  generator  adjusted  to  deliver  a  sinusoidal  8-period  voltage  pulse  with  amplitude  10 
V.  Since  the  function  generator  is  represented  by  a  50  fi  intrinsic  impedance,  a  25  resistor  was 
connected  in  series  with  it  to  obtain  a  total  source  impedance  of  75  Q.  In  the  same  maimer,  by 
connecting  a  25  SI  resistor  in  series  with  the  oscilloscope’s  50  fi  input,  a  receiver  load  resistance  of  75 
fl  was  obtained.  The  two  transducers  were  mounted  in  a  water  bath  of  dimensions:  1  m  x  0.5  m  x 
0.5  m  so  as  to  ensure  that  unintended  reflections  did  not  disturb  the  measurements.  Measurements 
as  well  as  numerical  results  were  obtained  corresponding  to  an  axial  transducer  distance  of  19  cm 
which  is  a  typical  transducer  distance  for  flow  meter  applications  in  the  flow  range:  6  -  40  m3/hour. 

The  dynamic  voltage  response  across  the  transmitter  electrodes  due  to  the  applied  generator 
voltage  in  Fig.  3a  is  given  in  Fig.  4a  (computed)  and  Fig.  4b  (measured).  Most  of  the  features 
seen  in  the  measured  response  are  also  found  in  the  computed  transmitter  voltage.  In  particular: 
a)  the  absolute  voltage  values  are  similar;  b)  the  first  period  in  the  response  is  characterized  by 
a  higher  amplitude  as  compared  to  the  subsequent  7  (i.e.,  N  —  1)  periods;  and  c)  as  soon  as  the 
generator  is  off,  the  transducer  continues  to  perform  damped  oscillations.  The  amplitude  of  the 
damped  oscillations  is,  however,  somewhat  higher  in  the  calculated  response.  This  is  due  to  mainly 
the  absence  of  mechanical  and  dielectric  losses.  Another  reason  is  that  the  transducer  is  a  three- 
dimensional  body  where  several  modes  arc  excited  (even  though  the  thickness  inode  is  by  far  the 
most  important  mode!)  and  stress  and  strain  cannot  exactly  be  represented  by  one  coordinate  only 
(z)  corresponding  to  the  quasi-plane  wave  assumption. 


(c) 


Figure  5:  Received  voltage  for  the  multilayer  structure  shown  in  Fig.  (2)  where  the  acoustic  coupling 
layer  thickness  is  chosen  to  be  zero,  a)  In  the  case  where  the  transmitter  is  ’short-circuited’  and  the 
receiver  is  open-circuited  ( Zi  =  0,Z2  =  oo);  b)  Loaded  by  75  Ohm  resistances  in  the  transmitter 
and  receiver  circuits  (Zi  =  75  Ohm,  Z2  -  75  Ohm);  c)  Measured  received  voltage  corresponding  to 
conditions  under  b). 


Computing  the  received  voltage  for  the  multilayer  transducer  structure  shown  in  Fig.  2  in  the 
case  where  no  series  resistance  appears  in  the  transmitter  circuit  (Zx  =  0)  and  at  the  same  time 
open-circuiting  the  receiver  transducer  ( Z2  =  oo),  i.e.,  in  the  limit  of  vanishing  electric  displacement 
between  the  receiver  electrodes,  results  in  the  received  voltage  pulse  depicted  in  Fig.  5a.  Since 
no  experimental  conditions  resemble  this  ideal  situation,  it  is  more  appropriate  to  include  a  series 
resistance  (which  is  chosen  to  be  a  75  Ohm  resistance  in  both  circuits;  a  typical  value  for  the 
electronics  in  flow  meter  applications)  when  comparing  theory  and  measurements.  In  doing  this,  the 
pulse  shown  in  Fig.  5b  is  computed.  It  is  evident  that  instead  of  two  consecutive  peaks  (Fig.  5a), 
one  peak  in  the  pulse  is  now  observed  (Fig.  5b),  and  this  latter  result  agrees  much  better  with 
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experimental  results  (Fig.  5c).  It  should  also  be  mentioned  that  the  computed  pulse  amplitude 
is  considerably  less  when  including  the  75  Ohm  resistances  (a  factor  of  4.77  is  found  between  the 
two  amplitudes).  In  order  to  have  a  common  standard  of  reference,  all  subsequent  results  arc 
normalized  with  respect  to  the  peak  value  of  Fig.  5b.  In  Fig.  5c,  measured  results  for  the  received 
voltage  response  is  given  corresponding  to  the  same  conditions  as  those  in  Fig.  5b.  The  general 
characteristics  are  the  same:  a)  the  peak  value  is  obtained  after  approximately  10  periods,  i.c., 
shortly  after  the  generator  is  in  off-mode;  and  b)  following  the  main  peak  of  the  pulse,  the  signal 
dies  away  characterized  by  a  time  constant  of  about  15  microsecs  (although  the  pulse  shows  somewhat 
stronger  damping  in  the  measured  case).  This  reflects  the  tendency  observed  for  the  voltage  across 
the  transmitter  electrodes  (Figs.  4a  and  4b). 

4  Conclusions 

A  complete  set  of  modeling  equations  for  one-dimensional  multilayer  ultrasound  transducers  is  pre¬ 
sented.  The  combination  of  Newton’s  Second  Law,  Poisson’s  equation,  and  the  definition  of  strain, 
allows  dynamic  equations  for  a  multilayer  transducer  to  be  determined.  Transducer  dynamics  be¬ 
come  completely  specified  by  imposing  continuity  of  normal  velocities  and  stress  across  material 
interfaces  as  well  as  the  electrical  condition  across  the  piezoccramic  electrodes.  Numerical  test  cases 
arc  discussed  for  a  transducer  consisting  of  a  piezoceramic  layer  and  a  half-wavelength  matching 
layer  of  stainless  steel.  Variations  in  transducer  reciprocal  signals,  accomodated  by  different  electri¬ 
cal  impedances  connected  to  the  transmitter  and  receiver,  are  finally  addressed,  and  good  agreement 
between  numerical  results  and  measurements  is  found. 
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1  Introduction 

Knowledge  of  acoustic  transmission  and  reflection  properties  of  construction  materials  has  a  great  importance 
for  the  silent  motion  of  vehicles.  The  purpose  of  this  study  was  to  experimentally  and  numerically  investigate 
an  Alberich  anechoic  rubber  coating  included  in  submarine  hulls.  Taking  favorable  acoustic  properties  into 
account  the  coating  can  be  manufactured  thinner  than  ordinary  coatings.  A  crucial  parameter  in  any  vehicle 
application  is  weight,  thus  by  using  the  thinner  anechoic  coating  the  weight  can  be  reduced. 

The  measurements  were  conducted  in  water  with  narrow  band  spherical  insonification  of  a  steel  plate,  coated 
with  the  anechoic  rubber.  Using  the  registered  data,  reflection  and  transmission  coeflicients  were  calculated 
based  on  integrated  direct,  reflected  and  transmitted  wave  packages.  The  considered  center  frequencies  were 
contained  in  the  kHz  region  and  the  used  experimental  setup  enabled  to  vary  the  angle  of  incidence.  Thus  the 
obtained  reflection  and  transmission  coefficients  could  be  expressed  as  function  of  center  frequency  and  angle 
of  incidence. 

A  computational  model  based  on  a  finite  element  approach  was  adopted  to  find  equivalent  elastic  parameters 
of  the  inhomogeneous  rubber  coating  in  terms  of  homogeneous  layers  by  matching  to  the  measured  data,  [1].  In 
the  model  the  plate  consists  of  multiple  homogeneous  layers  surrounded  by  water  halfepaces.  For  the  matching 
a  genetic  algorithm  was  employed  to  find  a  suitable  layer-parameter  combination,  [2). 

2  Measurements 

For  the  measurements  the  steel  plate  was  vertically  lowered  into  water.  The  plate  has  the  dimensions  2  by 
2  m  and  a  thickness  of  5  mm.  In  addition  the  plate  is  coated  with  a  5  mm  thick  Alberich  anechoic  rubber. 
The  rubber  contains  periodically  placed  cylindrical  air  inclusions  with  a  diameter  of  approximately  10  mm, 
resulting  in  the  anechoic  characterization.  For  reference  purposes  a  uncoated  steel  plate  was  used  for  verification 
of  the  method  and  the  equipment.  As  transmitters  and  receivers,  hydrophones  of  spherical  type  was  employed 
resulting  in  a  spherically  insonification  of  the  plate.  The  wave  package/source  signal  was  chosen  to  consist  of 
six  periods  of  a  cosine  function  with  Blackman  window  for  a  specified  center  frequency.  The  center  frequency 
ranged  from  10  kHz  to  30  kHz  with  a  step  of  5  kHz.  An  example  of  the  signal  is  displayed  in  Figure  1. 


Figure  1.  Left:  Transmitted  signal  with  the  Blackman  window. 

Right:  Measured  transmitted  signal  with  filtering.  Center  frequency  10  kHz. 

Furthermore  in  Figure  1  the  corresponding  measured  signal  is  displayed,  which  exhibits  a  slight  difference 
in  comparison  to  the  signal  specified  in  the  signal  generator.  This  small  perturbation  occurs  in  the  signal 
transition  between  hydrophone  and  water,  which  should  be  expected.  The  measured  signal  has  been  filtered 
and  then  used  as  input  for  the  numerical  computation. 
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Referring  to  Figure  2  the  experimental  setup  is  shown  with  scaled  distances.  The  hydrophone  where  the  wave 
package  is  transmitted  is  denoted  as  A.  With  hyrophones  placed  at  reflection  position,  B,  and  at  transmission 
position,  C,  resulting  in  a  angle  of  incidence  a,  time  signals  have  been  registered.  Using  the  direct  wave, 
the  transmitted  wave  and  the  reflected  wave  the  corresponding  reflection  and  transmission  curves  as  func¬ 
tion  of  center  frequency  and  angle  of  incidence  are  calculated.  The  angle  of  incidence,  a,  ranged  between 
8.13°  to  40.99°  degrees  and  was  varied  by  increasing/decreasing  the  depth  of  hydrophone  A  and  respectively 
decreasing/increasing  the  depth  of  hydrophones  B  and  C. 


Figure  2.  Experimental  setup  showing  orientations  of  the 
hydrophones  (A,B,C)  in  respect  to  the  investigated  plate. 

The  configuration  of  the  plate  and  the  hydrophones  has  been  chosen  with  care,  to  avoid  unwanted  interference 
from  later  arriving  field  components  in  the  registered  time  signals.  Such  components  are  for  instance  reflections 
from  the  water  surface  and  the  plate  edges. 

3  Computational  model 

In  the  computational  model  the  plate  is  composed  of  two  coupled  homogeneous  elastic  layers  surrounded  on 
both  sides  by  water  halfspaces.  A  cartesian  coordinate  system  (r ,z)  is  introduced  with  the  z-axis  orthogonal 
to  the  plate.  The  transmitter  is  represented  by  a  symmetric  point  source  on  the  z-axis,  exciting  a  cylindrically 
symmetric  wavefield  composed  of  compressional  waves  (P)  and  waves  with  polarization  parallel  to  the  z-axis 
(SV).  The  equations  of  motion  are  transformed  to  the  lateral  wavenumber  domain  (fc)  by  a  Hankel-Bessel 
transform. 

Ffm(h(fc))  =  f  h{r)Jm(kr)rdr  m  =  0,1,... 

Jo 

Here  dm  (At)  are  Bessel  functions  of  the  first  kind  and  k  is  the  wave  number  in  the  r  direction.  In  the  k  domain 
the  displacements  and  pressures  are  solutions  of  a  two  point  ordinary  differential  equation  boundary-value 
problem.  This  boundary-value  problem  is  solved  by  a  finite  element  approach  using  basis  functions  composed 
of  exact  local  solutions  in  the  layers.  The  finite  element  analysis  results  in  a  system  of  linear  equations. 

A  (fc)x  =  b(fc) 

The  equation  system  is  symmetric  and  block-tridiagonal  with  blocks  of  order  1  x  1  in  fluid  layers  and  2  x 
2  in  solid  layers.  The  corresponding  wavefield  components  in  the  (r,z)  domain  are  then  obtained  by  inverse 
Hankel-Bessel  transformation. 

-  J  H(kr)x(k)kdk 

The  contour  F  is  deviated  from  the  real  axis  to  an  equivalent  path  and  the  integral  is  evaluated  by  adaptive 
numerical  quadrature.  Transient  fields  are  obtained  by  Fourier  synthesis  of  their  mono-frequency  components. 
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A  comparison  in  the  time  domain  between  the  predicted  and  the  experimentally  observed  pulse  shape  in 
reflection  and  transmission  position  is  shown  in  Figure  3  for  the  uncoated  steel  plate.  The  parameters  of  steel 
and  the  water  halfspaces  are  assumed  to  be  known,  and  are  shown  in  Table  1.  In  Figure  3  data  is  normalized 
by  the  direct  pulse  arrival  at  the  reflection  position.  The  angle  of  incidence  is  24.3°  and  the  center  frequency 
is  20  kHz. 


Figure  3.  Upper-Left:  Reflection  position,  simulation.  Upper-Right:  Transmission  position,  simulation. 
Lower- Left:  Reflection  position,  measurement.  Lower-Right:  Transmission  position,  measurement. 

Note  that  both  curve  shapes  and  amplitudes  agree  for  the  uncoated  steel  plate  and  there  is  no  harmful  inter¬ 
ference  from  unwanted  field  components. 

4  Transmission  and  reflection  coefficients 

In  Figure  4  the  reflection  and  transmission  coefficients  for  the  uncoated  steel  plate,  with  parameters  in  Table 
1,  are  shown  as  function  of  incidence  angle  at  a  center  frequency  of  20  kHz. 


ANGLE  ANGLE 

Figure  4.  Left:  Reflection  coefficient.  Right:  Transmission  coefficient. 

The  average  deviations  for  the  reflection  and  transmission  coefficients  are  ±0.21  dB  respectively  ±0.42  dB, 
which  is  satisfactory. 

Turning  to  the  coated  steel  plate  the  elastic  parameters  of  the  homogeneous  rubber  layer  should  be  matched 
to  the  measured  data.  For  the  matching  a  genetic  algorithm  is  employed.  The  algorithm  finds  a  suitable 
parameter  combination  for  the  rubber  layer  in  the  parameter  space  within  specified  boundaries.  In  this  case 
the  parameter  space  consists  of  the  density,  the  P  wave  velocity,  the  S  wave  velocity  and  the  thickness  for 
the  layer.  The  parameter  sets,  individuals,  are  represented  by  bit-strings,  ’chromosomes’,  unique  for  each 
individual.  Initially  a  finite  number  of  individuals  are  chosen  at  random,  constituting  the  first  generation. 
Then,  a  simulated  evolution  proceeds  by  applying  mating,  reproduction  and  mutation  operators  iteratively. 
Eventually  a  best  fitted  individual  is  produced  for  a  specific  number  of  iteration  steps. 

For  the  matching  the  fitness  function  has  been  chosen  to  be  the  the  sum  of  the  average  difference  between 
the  measured  and  the  predicted  reflection  and  transmission  coefficients,  in  dB  scale.  The  measured  coefficients 
are  based  on  maximum  amplitudes  of  the  the  observed  wave  packages.  For  the  predicted  coefficients  a  simplifi¬ 
cation  has  been  introduced  by  using  computed  mono-frequency  fields,  thus  making  the  computations  less  time 
consuming. 
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The  resulting  material  parameters  axe  shown  in  Table  1.  Describing  the  inhomogeneous  coating  by  one  homo¬ 
geneous  layer  resulted  in  a  thickness  of  50.2  mm,  consequently  an  approximately  10  times  greater  thickness. 


Material 

P-vel  m/s 

P-atten  db/A 

S-vel  m/s 

S-atten  db/A 

Water 

1000.0 

1@H  jgili 

- 

- 

- 

Steel 

7700.0 

- 

3230.0 

_ 

Rubber 

544.6 

330.2 

1.0 

84.9 

50.0 

Table  1.  Material  parameters  of  the  surrounding  medium  and  the  plate. 


In  figure  5  the  corresponding  reflection  and  transmission  coefficients  for  the  coated  steel  plate  is  displayed  at 
a  center  frequency  of  20  kHz. 


iao  iao  zan  zs.o  30.0  35.0  4a  o  45.0 

ANGLE 
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Figure  5.  Left:  Reflection  coefficient.  Right:  Transmission  coefficient. 


The  average  deviations  for  the  reflection  and  transmission  coefficients  axe  ±0.27  dB  respectively  ±0.48  dB, 
which  at  this  stage  is  satisfactory. 


5  Conclusions  and  future  work 

For  the  uncoated  steel  plate  good  agreement  between  the  observed  and  the  theoretically  predicted  pulse  shapes 
has  been  achieved,  see  Figure  4.  For  the  coated  steel  plate  the  agreement  is  satisfactory,  see  Figure  5,  but 
the  thickness  of  the  homogeneous  rubber  layer  is  approximately  10  times  greater  than  the  real  thickness 
of  the  coating.  The  future  objective  is  to  obtain  the  simplest  equivalent  material  for  the  coating  with  the 
prerequisite  to  reduce  the  overall  thickness  of  the  composite  layer.  For  this  purpose  several  layers  describing 
the  inhomogeneous  rubber  could  be  introduced.  Matching  all  frequencies  simultaneously  will  be  of  particular 
interest,  since  the  resulting  parameter  combination  could  be  used  as  a  model  of  the  rubber  layer  in  a  broad 
frequency  range. 
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Abstract 

High-resolution  sonar  at  long  range  is  desirable  for  future  mine  hunting  and  classification.  Synthetic  Aperture 
Radar  (SAR)  is  now  used  routinely  to  achieve  high  resolution  at  long  distances,  so  it  is  natural  to  look  to  Synthetic 
Aperture  Sonar  (SAS)  for  a  solution  to  the  problem.  Unfortunately  the  operational  application  ofSAS  has  been  very 
limited.  The  low  velocity  of  sound  compared  with  radar  signals  creates  a  practical  problem.  Variability  of  the 
propagating  medium  and  uncertainty  in  position  of  the  sonar  platform  create  problems  in  achieving  the  necessary 
phase  coherence  for  aperture  synthesis.  SAS  image  resolution  depends  on  how  well  these  effects  can  be  compen¬ 
sated.  In  SAR,  correction  of  platform  position  from  the  radar  data  itself  is  termed  "autofocus".  The  term 
"autopositioning”  is  preferred  for  SAS,  since  image  focus  is  seldom  used  directly.  The  paper  describes  an  experi¬ 
ment  using  the  Remotely  Operated  Vehicle  (ROV),  PLUMS,  and  a  wideband  sonar  with  a  multiaperture  receiver 
array.  Artificial  targets  were  placed  on  the  seabed  about  100  m  from  the  ROV  track.  A  beacon  positioning  system 
with  two  transmitters  on  the  seabed  was  used  for  initial  positioning.  The  availability  of  a  physical  aperture  image 
from  each  transmitted  pulse  greatly  simplified  the  autopositioning  problem.  DPC  (Displaced  Phase  Centre)  meth¬ 
ods  were  investigated,  based  on  both  seabed  and  target  image  correlation.  The  simplest  method  investigated  wm  a 
variant  of  Prominent  Point  Processing  as  used  in  SAR.  The  paper  provides  some  theoretical  background  to  these 
different  methods,  together  with  corrected  images  of  the  experimental  targets. 


1 1ntroduction 

The  development  of  new  mines,  which  can  move  and 
hit  ships  at  some  distances,  has  made  mine  hunting  more 
risky  for  the  Mine  Counter  Measure  (MCM)  ships.  The 
capability  to  detect  and  classify  mines  at  a  distance  is 
therefore  highly  desirable,  especially  in  non-war  situa¬ 
tions  where  personnel  losses  cannot  be  tolerated.  There 
are  at  least  two  ways  to  achieve  this  capability.  One  is  to 
mount  the  sonar  on  a  Remotely  Operated  Vehicle  (ROV), 
which  can  then  be  moved  close  to  the  mine  without  en¬ 
dangering  the  MCM  ship.  Another  way  is  to  improve 
the  spatial  resolution  of  the  sonar,  so  that  classification 
can  be  done  at  a  safe  distance.  Higher  resolution  requires 
a  bigger  sonar  aperture,  but  there  are  practical  and  eco¬ 
nomical  limits  for  the  physical  size  of  an  MCM  sonar, 
and  current  designs  are  coming  close  to  this.  There  is 
therefore  considerable  interest  in  the  use  of  Synthetic 
Aperture  Sonar  (SAS)  to  increase  the  effective  aperture 
size  for  classification  purposes.  Range  resolution  is  also 
important  for  SAS  imaging,  partly  for  the  final  image 
quality,  and  partly  for  autopositioning  along  the  way. 
We  have  therefore  investigated  the  use  of  SAS  for  clas¬ 
sification  of  MLOs  (mine-like  objects)  at  100m  range, 
in  an  experiment  using  an  ROV-mounted  wideband  so¬ 


nar.  The  wideband  SAS  processing  method  consists  of 
combining  successive  physical  aperture  images  as  de¬ 
scribed  in  [I],  This  is  a  conceptually  simple  time-do¬ 
main  method  for  wideband  signals,  avoiding  the  band¬ 
width  and  other  limitations  of  many  FFT-based  process¬ 
ing  techniques.  The  problem  of  slow  processing  in  the 
time  domain  will  be  discussed  later  in  the  paper. 

2  The  Field  Experiment 

The  experiment  was  performed  at  the  FOI  test  site 
Djupviken  in  October  1999.  Djupviken  is  a  semi-en¬ 
closed  bay  in  the  southern  Stockholm  Archipelago, 
where  the  seabed  consists  of  soft  postglacial  clay  and 
mud.  At  this  site,  there  is  a  pontoon  laboratory,  which  is 
completely  floating  in  10  m  of  water.  The  seafloor  slopes 
downwards,  from  a  depth  of  10m  where  the  ROV  was 
operated,  to  20m  depth  at  200m  range  out  in  the  bay. 
The  seabed  is  fairly  smooth  and  free  from  vegetation. 
Two  transmitters,  “beacons”,  were  placed  on  the  seabed 
and  used  to  determine  the  position  of  the  platform  by 
triangulation.  The  beacons  were  placed  80m  from  the 
ROV  track,  about  40m  apart  (Figure  1). 
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Figure  1.  The  objects  at  the  experimental  site. 


Three  targets  were  positioned  about  100m  from  the 
ROV  track,  and  the  target  area  was  delimited  by  four 
comer  reflectors.  The  targets  were  two  MLOs,  the  “Small 
Target”  with  dimensions  0.4x0. 4x1. 6m  and  the  “Large 
Target”  (07x0.7x1 ,8m),  together  with  the  “T-Target”,  de¬ 
signed  to  assess  spatial  resolution  in  the  SAS  image.  This 
consisted  of  a  T-formed  metal  frame  with  plastic  foam 
balls  attached  on  short  rods  above  it  (Figure  2).  During 
processing  another  target  was  discovered  at  closer  range, 
which  turned  out  to  be  a  concrete  slab  sunk  in  the  seabed. 

The  sonar  was  mounted  on  the  PLUMS  (Platform 
for  Underwater  Measurement  Systems)  ROV,  (Figure  3), 
designed  by  the  Swedish  Defence  Research  Agency. 
PLUMS  is  an  open  stainless  steel  pipe  frame  within 
which  different  types  of  equipment  can  be  mounted. 
Heading  and  attitude  are  measured  using  a  three-axis  rate 
gyro  aided  by  a  magnetic  compass  and  two  pendulums. 
Roll,  pitch  and  yaw  angles  can  be  controlled  to  errors  of 
less  than  one  degree.  Depth  can  be  kept  constant  to  within 
0.1m  down  to  the  maximum  depth  of  70  m. 


Figure  3.  The  ROV  PLUMS. 


Two  different  transmitters  were  used  for  imaging, 
an  EDO  side-scan  array  and  a  200  kHz  single  trans¬ 
ducer  transmitter,  here  termed  Tx200.  The  receiver  ar¬ 
ray  consisted  of  two  identical  1 6-element  sub-arrays  with 
a  spacing  of  15mm.  For  the  Djupviken  experiments,  they 
were  mounted  side  by  side  to  form  a  32-element  Uni¬ 
form  Linear  Array  (ULA).  The  DAIM  Receiver  has  an 
operational  bandwidth  from  50  kHz  to  100  kHz,  which 
is  extremely  high.  It  is  difficult  and  costly  to  achieve 
uniform  receiver  response  over  a  large  sector  in  azi¬ 
muth  and  a  wide  bandwidth.  However  in  the  analysis, 
the  received  signals  are  always  first  matched  filtered 
with  a  computed  kernel.  The  kernel  is  based  on  a  recor¬ 
ded  replica  of  the  source  signature  and  its  Hilbert  trans¬ 
form  and  can  be  adjusted  to  give  equal  response  for  all 
receiver  channels  [2],  Software  calibration  reduces  the 
need  for  uniform  response  from  the  receiver  hardware. 

In  our  experiments,  the  beacons  and  the  EDO  array 
used  a  linear  chirp  from  60kHz  to  120kHz,  while  the 
Tx200  chirped  from  120kHz  to  240kHz.  Unfortunately 
the  receiver  response  was  not  constant  within  this  fre¬ 
quency  span,  resulting  in  a  usable  bandwidth  of  only 
about  twenty  percent  of  the  nominal  bandwidth.  Match 
filtering  also  decreases  the  effective  bandwidth,  since 
the  spectrum  is  in  some  sense  squared. 


3  Beacon  Positioning 

3.1  Introduction 

In  SAS,  image  resolution  depends  on  positioning  accu¬ 
racy.  Along-track  errors  of  the  order  of  the  acoustic 
wavelength  and  even  smaller  across-track  errors  lead  to 
evident  phase  incoherence,  giving  blurred  focus  and  high 
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(3.1) 


sidelobe  levels.  Initial  positioning  estimates  can  be  ob¬ 
tained  by  dead  reckoning  with  a  high-quality  platform. 
However  the  PLUMS  vehicle  was  not  equipped  for  pre¬ 
cision  inertial  navigation,  so  beacon  positioning  was  used 
instead.  Two  wideband  transmitters  mounted  on  the 
seabed  were  used  as  baseline.  In  previous  rail-based 
experiments  using  the  DAIM  receiver  array,  very  accu¬ 
rate  range  estimates  were  made  using  a  shore-mounted 
wideband  transmitter,  and  it  was  hoped  to  repeat  this 
experience.  Look-angles  up  to  1  radian  are  possible  with 
the  DAIM  receiver,  allowing  positioning  along  an  ex¬ 
tended  track.  However  beacon  positioning  is  an  opera¬ 
tional  inconvenience,  and  inertial  platforms  are  expen¬ 
sive,  so  a  main  objective  of  the  experiment  was  to  in¬ 
vestigate  autopositioning  methods,  where  the  platform 
position  is  corrected  from  the  acoustic  data  itself. 

3.2  Position  Estimation 

The  precise  separation  between  the  beacons  was  not 
known.  Hence  the  positioning  computation  was  per¬ 
formed  in  three  stages.  First  range  and  look-angle  was 
estimated  for  all  pings  to  one  beacon.  Range  and  look- 
angle  was  then  estimated  for  all  pings  to  the  second  bea¬ 
con.  The  look-angles  as  seen  from  the  receiver  array 
were  estimated  using  wideband  MUSIC  [3,4],  and  ena¬ 
bled  beacon  separation  to  be  estimated  by  the  cosine 
law.  Vehicle  position  could  then  be  determined  by  trian¬ 
gulation  using  the  two  ranges  and  beacon  separation. 
The  angular  information  also  allowed  vehicle  heading 
to  be  determined  with  more  precision  than  was  given  by 
the  on-board  instrumentation. 

The  matched  filtered  signal  was  first  windowed  to 
remove  any  surface  or  bottom  reflections,  which  also 
eliminated  signals  from  the  second  beacon  and  other 
sources.  The  platform  moved  at  a  slow  steady  speed  so 
the  window  was  chosen  around  the  range  determined 
from  the  previous  ping.  The  MUSIC  look-angle  allowed 
beamforming  in  the  beacon  direction  followed  by  range 
estimation.  Signal  travel  times  were  determined  either 
from  the  peak  amplitude  of  the  compressed  signal,  or 
by  correlation  with  a  matched  filtered  replica  of  the  trans¬ 
mitted  pulse. 

3.2  Navigation 

Range  and  angle  estimates  are  made  in  the  plane  con¬ 
taining  the  receiver  and  beacon.  Triangulalion  is  required 
in  a  horizontal  plane,  so  these  ranges  and  angles  are  first 
resolved  into  the  horizontal  plane.  The  baseline  between 
the  two  beacons  can  then  be  computed  from  the  cosine 
law, 


D2  =  r2  +  r22  -  2 rxr2  cos(0j  +  02 ) 

where  r  and  theta  are  the  horizontal  range  and  look- 
angles  to  each  beacon  (Figure  4).  Separation  between 
the  beacons  is  estimated  as  the  mean  separation  over  all 
pings.  Platform  position  is  determined  in  the  beacon 
coordinate  system,  with  origin  midway  between  the 
beacons.  The  x-axis  points  along  the  baseline  while  the 
y-axis  is  orthogonal  to  the  baseline  pointing  towards 
the  vehicle  (Figure  5).  The  platform  position  is  given  as 


and 


(3.2) 


(3.3) 


Vehicle  heading  is  determined  from  the  bearing  and 
look-angle  to  both  beacons,  and  then  averaged. 

In  the  Djupviken  experiment,  the  beacon  baseline 
was  quite  short  compared  with  the  range  to  the  vehicle, 
while  the  vehicle  path  was  almost  parallel  to  the  base¬ 
line.  The  ensuing  along-track  navigation  errors  were 
greater  than  the  cross-track  errors.  Estimates  made  from 
the  imagery  suggest  an  along-track  random  position  er¬ 
ror  of  74  mm,  and  an  across-track  random  error  of  1 8mm. 
These  are  not  absolute  errors  but  the  errors  relevant  to 
phase  coherence  along  the  aperture. 

Beacon  2 


r2 


Platform 


Figure  4.  The  platform 
and  the  two  beacons 
projected  on  a  horizon¬ 
tal  plane. 


Figure  5.  The  platform, 
the  beacons  and  the 
beacon  coordinate 
system. 
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4  Autopositioning 

4.1  Introduction 


In  the  SAR  community,  and  in  much  of  the  SAS  litera¬ 
ture,  aulopositioning/aulofocus  is  entirely  concerned 
with  cross-track  navigation  error.  This  is  because  phase 
coherence  is  so  much  more  sensitive  to  cross-track  er¬ 
ror.  Here  the  along-track  errors  were  too  large  to  ignore, 
so  methods  which  could  correct  for  both  components  of 
error  were  investigated.  A  review  of  SAR  aulofocus 
methods  is  contained  in  [5],  Phase  gradient  methods  can 
detect  track  curvature,  but  are  less  helpful  for  the  ran¬ 
dom  ping-to-ping  platform  errors  which  are  a  feature  of 
SAS.  However  one  method,  Prominent  Point  Processsing 
can  be  implemented  in  a  particularly  simple  way  using 
a  multi-channel  receiver  array.  This  will  be  described  in 
Section  4.2.  Displaced  Phase  Centre  (DPC)  methods 
appear  to  have  been  developed  specifically  for  sonar 
[6],  Leading  and  trailing  subarrays  are  selected  from  the 
full  receiver  array  in  such  a  way  that  the  phase-centre, 
defined  as  the  midpoint  between  transmitter  and  receiver, 
remains  stationary  between  consecutive  pings.  In 
SherrilFs  implementation,  seabed  echoes  from  consecu¬ 
tive  pings  are  correlated,  and  the  phase  displacement 
gives  cross-track  (strictly  cross-heading)  error.  In  our 
method,  DPC  images  of  the  seabed  are  correlated  in¬ 
stead  of  DPC  echoes.  Section  4.3  shows  that  this  method 
is  capable  of  giving  both  along-  and  cross-heading  er¬ 
ror.  A  similar  method  can  be  used  with  target  images. 

SAS  processing  is  usually  carried  out  in  the  FFT 
Domain.  However  when  a  long,  wideband,  receiver  ar¬ 
ray  is  used,  and  random  autopositioning  corrections  are 
required,  it  is  more  convenient  to  generate  the  SAS  im¬ 
age  by  summing  physical  aperture  images  derived  from 
successive  pings,  leading  to  a  lime  domain  back-pro¬ 
jection  approach.  Autopositioning  is  then  also  based  on 
physical  aperture  images  generated  at  each  platform 
position.  The  disadvantage  with  this  method  is  slow 
processing,  but  faster  time-domain  algorithms  are  now 
becoming  available  [7], 

4.2  Prominent  Point  Processing 

This  is  conceptually  the  simplest  method.  As  noted  ear¬ 
lier,  SAS  imaging  is  carried  out  by  adding  physical  ap¬ 
erture  images,  generated  at  each  point  along  the  aper¬ 
ture.  If  a  strong  isolated  reflector  is  present  in  the  target, 
this  gives  rise  to  a  point-spread  function  in  each  image 
with  peak  brightness  at  the  pixel  corresponding  to  the 
reflector  location.  Phase  coherence  for  this  reflector  is 
then  forced  by  phase-shifting  each  physical  aperture 
image  so  that  the  phase  at  the  selected  pixel  location  is 


zero.  This  phase  shift  corrects  for  both  along  and  across- 
heading  errors  resolved  along  the  line  of  sight  to  the 
selected  pixel  in  the  target  area,  and  thus  for  the  whole 
image  if  the  area  is  sufficiently  small.  The  same  pixel 
location  must  be  used  throughout,  since  the  brightest 
pixel  may  move  from  one  ping  to  the  next.  The  strategy 
used  with  recorded  data  was  to  select  the  location  of  the 
brightest  pixel  in  an  image  generated  from  the  middle 
of  the  synthetic  aperture. 

Only  one  isolated  strong  reflector  should  be  present 
at  the  selected  location.  Otherwise  phase  interference 
between  neighbouring  reflectors  causes  errors.  Also  the 
uncorrecled  navigation  must  be  sufficiently  good  to  se¬ 
lect  the  same  point  in  the  image  each  time.  However  the 
image-based  method  requires  no  interaction,  and  a  sin¬ 
gle  phase-shift  at  mid-ftequency  could  be  used  to  han¬ 
dle  the  wideband  echoes  [8]. 

4.3  Seabed  Image  Correlation 

The  idea  is  to  correlate  images  of  the  same  area  of  seabed 
from  consecutive  pings,  generated  using  DPC  subarrays 
chosen  from  the  full  receiver  array.  If  there  is  no  dis¬ 
placement,  then  the  consecutive  images  will  correlate 
perfectly.  However  if  the  phase-centre  is  displaced  by  a 
small  amount  on  the  seabed,  then  complex  images  of 
seabed  or  target  are  phase-shilled  by  an  angle  propor¬ 
tional  to  platform  displacement  [9], 

For  seabed  autopositioning,  the  coordinates  of  the 
image  frame  move  with  the  phase-centres  of  the 
subarrays,  as  given  by  the  navigation  data.  The  image 
area  should  lie  directly  abeam  of  the  displaced  phase 
centres  for  consecutive  pings,  designed  to  be  almost 
coincident  in  the  along-track  direction.  Consider  the 
following  simplified  model.  The  sonar  platform  travels 
with  constant  heading  along  the  x-axis,  with  the  receiver 
axis  pointing  in  the  y-direclion.  The  seabed  is  defined 
by  a  random  distribution  of  reflectors  on  a  plane,  hori¬ 
zontal,  surface.  For  simplicity,  this  is  assumed  to  be 
coplanar  with  the  vehicle,  an  assumption  not  made  in 
actual  processing.  Define  an  H  x  K  rectangular  image 
frame  on  this  surface  centred  at  [X,Y],  Suppose  that  a 
ping  occurs  when  the  phase  centre  of  the  array  lies  at 
[0,0],  A  physical  aperture  complex  image  is  made  of  the 
reflectors  visible  within  the  frame,  with  coordinates  [x\ 
y’]  centred  on  the  middle  of  the  frame.  Suppose  Y  is 
sufficiently  large  for  all  points  in  the  image  to  lie  in  the 
far  field  of  the  array,  and  that  Y»X,  Y»H,  Y»K.  Then 
all  angles  a  are  small  enough  for  sin  a=a.  Consider  a 
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reflector  located  at  x’=  a,  y’=b  within  the  frame.  Then 
the  psf  is  a  circular  arc  through  [a,b],  approximating  a 
straight  line,  making  angle  a  =  -atan((X+a)/(Y+b))  with 
the  X-axis,  ie 


X  +  a 
a  =  — — — 

Y 


(41) 


A  second  ping  is  transmitted  when  the  phase-centre  has 
moved  to  [u,v].  There  is  no  navigation  error,  so  the  im¬ 
age  frame  also  moves  by  [u,v].  The  reflector  moves  to 
[a-u,  b-v]  within  the  frame.  The  change  in  a  is  negligi¬ 
ble.  The  two  images  are  now  correlated  in  the  column, 
or  Y,  direction.  Suppose  that  one  of  the  columns  coin¬ 
cides  precisely  with  x’  =  a.  The  psf  in  the  second  image, 
centred  on  [a-u,  b-v]  intercepts  this  column  at 


u(X  +  a) 

y’=b  -v-uima  ~b  -v - — - 


(4.2) 


Then  image  correlation  for  this  column  will  show  a 
displacement  given  by 


Ay’(o) 


u(X  +  a) 
7 


(4.3) 


If  there  are  many  reflectors  in  the  scene,  then  the 
displacements  for  all  the  columns  across  the  image  lie 
on  a  straight  line  given  by 


Ay’(*’) 


X  +  x’ 
Y 


(4.4) 


The  slope  of  this  correlation  line 


(4.5a) 


depends  on  the  along-heading  displacement,  u.  Write 
Q=XfY~  (minus  look-angle)  for  X«Y.  The  intercept  at 
x’=0  is  given  by 


Ay’(o)  =  -v-0w  (4.5b) 

Equations  (4.5)  can  be  inverted  to  obtain  [u,v] ,  the  true 
movement  of  the  phase-centre  from  the  slope  and  inter¬ 
cept  of  the  image  correlogram.  The  navigation  error  [g,h] 
is  the  difference  between  the  [u,v]  and  the  displacement 
given  by  the  navigation  file.  To  obtain  the  integrated 
track  error,  these  ping-to-ping  errors  [g,h]  must  be 
summed  over  the  pings  along  the  synthetic  aperture.  The 


above  analysis  was  confirmed  by  simulation  studies  us¬ 
ing  random  distributions  of  point  reflectors  to  represent 
the  seabed. 

The  above  arithmetic  does  not  rely  on  stabilizing  the 
phase  centre  against  platform  movement,  but  DPC  is 
needed  to  retain  cross-correlation  between  images  when 
the  platform  moves  a  substantial  distance  between  pings. 
In  our  experiments,  the  correlation  levels  between  seabed 
images  was  poor,  so  a  lot  of  work  went  into  extracting 
the  line  fit  required  for  Equation  4.4.  Images  are  already 
broken  into  columns  by  the  frame  definition.  They  w  ere 
further  divided  into  overlapping  sections  in  the  row  di¬ 
rection.  Each  section  of  each  column  was  then  corre¬ 
lated  separately  with  the  matching  section  in  the  partner 
image.  Correlation  of  each  pair  of  sections  yields  four 
quantities  -  peak  absolute  ccf,  ratio  of  peak  ccf  to  next 
highest  peak,  estimated  Y  displacement,  cross-power. 
The  correlation  is  rejected  if  peak  ccf  or  peak  ratio  are 
too  low  or  Y  displacement  is  too  high.  The  result  of  im¬ 
age  correlation  is  a  pair  of  matrices  of  Y-displacements 
and  cross-powers,  with  rejected  sections  indicated  by 
NaN  entries.  A  Hough  Transform,  weighted  by  cross¬ 
power,  was  then  be  used  to  find  the  best  line  fit. 

4.4  Target  Image  Correlation 

A  prominent  landmark  such  as  a  rock  outcrop  can  be 
used  in  a  similar  manner  to  the  seabed,  but  leaving  the 
image  frame  fixed  as  the  platform  moves.  The  arithme¬ 
tic  is  slightly  different  but  the  principle  is  the  same.  The 
prominent  object  can  be  the  target  itself.  However  the 
target  is  too  small  to  be  resolved  in  the  physical  aper¬ 
ture  image  (otherwise  SAS  would  not  be  necessary),  so 
will  not  appear  as  a  random  reflector  distribution.  Us¬ 
ing  a  point-target  model,  the  psf  is  a  sector  of  a  circular 
arc  centred  on  the  phase-centre  at  each  platform  posi¬ 
tion.  A  little  arithmetic  leads  to  the  expected  result  that 
the  y-inlercept  in  (4.5b)  gives  a  combined  correction  for 
navigation  error  resolved  in  the  direction  of  the  target. 
This  correction  was  used  for  the  results  shown  in  the 
paper.  Subject  to  certain  conditions,  it  is  also  possible 
to  separate  out  along  and  across  heading  error  compo¬ 
nents,  but  this  was  not  done.  The  theory  was  confirmed 
by  simulation  studies,  including  the  simulated  T-Target 
showing  that  the  point  target  model  was  appropriate  in 
this  case. 
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5  SAS  Processing 

Prominent  Point  Processing 

In  spite  of  its  theoretical  limitations,  this  method  worked 
consistently  well  with  all  targets  imaged.  Examples  of 
PPP  compensation  are  shown  in  Figure  6  (T-target)  and 
Figure  1 2  (Large  target).  Of  course  any  detected  target  is 
likely  to  be  a  prominent  object  in  the  scene,  but  it  is  less 
obvious  that  the  target  will  contain  prominent  points, 
which  are  satisfactory  for  the  algorithm  at  the  resolution 
required  for  classification. 

Seabed  Autopositioning 

Conditions  for  seabed  autopositioning  were  difficult  at 
Djupviken.  Due  to  the  sloping  seabed  no  seabed  echoes 
were  received  from  target  range  using  the  200  kHz  trans¬ 
mitter,  and  no  seabed  echoes  at  all  using  the  EDO  array, 
which  was  unfortunate  since  the  T-Target  was  only  vis¬ 
ible  with  that  transmitter.  At  200  kHz  useful  seabed  ech¬ 
oes  were  received  out  to  a  range  of  35  m  from  the  PLUMS 
track,  but  there  were  severe  problems  with  surface  mul¬ 
tiples,  since  the  DAIM  receiver  has  a  wide  elevation 
beam  pattern.  There  were  also  problems  with  periodi¬ 
cally  corrupted  echoes,  which  made  it  difficult  to  find  a 
long  enough  sequence  of  pings  for  SAS  imaging. 

In  the  absence  of  true  position,  one  test  of  seabed 
autopositioning  is  to  compare  position  estimates  using 
images  at  different  ranges.  Comparative  results  at  200 
kHz  for  mid-ranges  21m  and  25m  are  shown  in  Figure 
7.  The  results  for  mid-ranges  21m  and  28m  are  shown  in 
Figure  8.  The  comparison  looks  quite  encouraging  apart 


Figure  6.  SAS  image  of  T-target  with  EDO  array  and  PPP 
compensation. 


from  a  few  spikes.  However  when  the  summed  ping-to¬ 
ping  errors  are  compared  in  Figure  9,  it  can  be  seen  that 
the  corrected  tracks  jump  apart,  although  they  remain 
parallel  for  short  sequences  at  a  time.  This  shows  the 
susceptibility  of  incremental  methods  to  occasional  large 
errors.  The  result  of  Seabed  Autopositioning  is  illustrated 
by  Figures  10-1 1  where  the  “Large  target”  is  imaged  with 
and  without  Seabed  Autopositioning. 


Figure  7.  Seabed  Autopositioning  at  200  kHz.  Ping-to-ping 
cross-track  error  comparison  for  mid-ranges  21  m  (solid 
line)  and  25  m  (dashed  line). 


Figure  8.  Seabed  Autopositioning  at  200  kHz.  Ping-to-ping 
cross-track  error  comparison  for  mid-ranges  21  m  (solid 
line)  and  28  m  (dashed  line). 


Figure  9.  Seabed  Autopositioning  at  200  kHz.  Summed 
ping-to-ping  cross-track  error  comparison  for  mid-ranges 
21  m  (lower  line)  and  25  m  (upper  line). 
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Figure  10.  SAS  image  of  "Large  target”  withTx200,  Seabed 
Autopositioning  and  PPP  compensation. 


Figure  11.  SAS  dB  image  of  "Large  target”  with  Tx200 
and  PPP  compensation. 


Target  Autopositioning 

Good  image  correlations  were  obtained  at  both  frequen¬ 
cies,  using  the  method  outlined  in  Section  4.3.  The  re¬ 
sult  for  the  T-Target  was  quite  similar  to  that  of  the  PPP 
method. 

Comparison 

Figures  11  and  12  show  images  of  the  Large  Target  us¬ 
ing  each  transmitter.  Theoretically  the  Tx200  should  pro¬ 
duce  better  resolution  but  to  achieve  the  higher  resolu¬ 
tion  positioning  must  be  good  enough.  Since  the  posi¬ 
tioning  accuracy  is  equivalent  in  both  cases  the  best  im¬ 
age  is  actually  produced  with  the  EDO  array. 


Figure  12.  SAS  dB  image  of  "Large  target"  with  EDO  and 
PPP  compensation. 


6  Conclusions 

The  beacon  positioning  system  using  the  wideband  re¬ 
ceiver  array  worked  well,  but  accuracy  was  inadequate 
for  synthetic  aperture  imaging.  However  after 
autoposition  correction,  sufficiently  good  SAS  images 
for  target  classification  were  obtained. 

In  spite  of  theoretical  limitations.  Prominent  Point 
Positioning  worked  consistently  well.  This  algorithm 
proved  particularly  simple  to  implement  using  the  multi¬ 
channel  receiver  array.  An  alternative  method  based  on 
correlating  DPC  physical  aperture  images  of  the  target 
also  worked  well.  This  is  important,  because  the  DPC 
method  is  more  general.  Autoposition  by  correlation  of 
DPC  seabed  images  has  been  shown  theoretically  and 
by  simulation  to  yield  both  along  and  across-heading 
errors,  but  experimental  results  were  disappointing.  Fur¬ 
ther  work  is  needed  to  decide  whether  these  reflect  a 
fundamental  sensitivity  of  the  method  to  the  seabed  en¬ 
vironment. 

For  future  research,  there  are  a  number  of  practical 
steps  which  can  be  taken  to  improve  the  experimental 
equipment.  One  goal  for  this  research  is  to  develop  ro¬ 
bust  autopositioning  methods,  which  are  adequate  for 
SAS  imaging  in  combination  with  simple  dead-reckon¬ 
ing  navigation  on  the  platform. 
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Simulation  of  nonlinear  fields  from  2D  ultrasound  ^ 
transducers  and  a  new  secondary  grating  lobe 

phenomenon 

Helge  Fjellestad  and  Sverre  Holm 


Abstract— la  this  study  we  hare  used  Burgers’  equation  to  study  the  field 
from  sparse  ultrasound  arrays  where  the  full  array  had  50  by  50  elements. 
The  two-way  sidelobe  level  was  reduced  by  20-25  dB  compard  to  the  fun¬ 
damental  mode.  We  have  also  found  that  a  new  grating  lobe  appears  at  the 
second  harmonics  at  half  the  sine  of  the  angle  of  the  fundamental  grating 
lobe. 


I.  Introduction 

In  this  study  we  wanted  to  simulate  the  two-way  response 
from  sparse  arrays  when  imaging  in  a  nonliner  medium.  The 
purpose  was  to  compare  the  performance  under  nonlinarity  with 
that  in  a  linear  medium  as  found  in  [1], 

II.  Theory 

We  have  used  Burgers’  equation  [2]  to  describe  the  acoustic 
field: 


du  _  (3uj0  du  d2u 

dz  <?Q  U  dr  dr 2 


(1) 


In  this  expression,  c0  is  the  nominal  speed  of  sound  in  the  me¬ 
dium  (speed  of  infinitesimal  bulk  waves),  /?  =  1  +  ^  where 
§  is  the  ratio  of  the  first  two  terms  in  the  nonlinear  pressure 
-  density  relation  for  the  medium,  and,  finally,  T  is  a  constant 
related  to  the  thermo-viscous  dissipation  of  the  medium.  Also, 
t  =  tut  —  kz. 

The  expression  can  be  modified  for  numerical  implementa¬ 
tion  [3]: 


17  =  „  +  f)  (2) 

^  °  m=l  m=n 

The  infinite  series  has  to  be  truncated  to  a  finite  number  of 
harmonics,  N.  The  simulator  is  based  on  the  principles  of  [4] 
and  with  modifications  from  [3]  in  order  to  generalize  it  from  a 
circular  source  geometry  to  a  rectangular  geometry.  In  the  sim¬ 
ulator,  nonlinear  propagation  is  taken  care  of  through  the  trun¬ 
cated  eq.  2  and  diffraction  is  handled  from  one  depth  to  another 
through  the  angular  spectrum  method  [5].  Attenuation  can  also 
be  introduced  in  the  diffraction  step,  although  not  used  in  our 
study. 
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Strandveien  4,  P.O.Box  354,  N-1326  Lysaker,  Norway.  Sverre  Holm  is  with  the 
Department  of  Informatics,  University  of  Oslo,  P.O.Box  1080  Blindern,  N-0316 
Oslo,  Norway. 


Transducer: 
Beam  type 
Radius 
Initial  pressure 
Frequency 
Elevation  focus 
Azimuth  focus 
Wavelength 

CW 

7.5  mm 

600  kPa 

3  MHz 

80  mm 

80  mm 

0.5133  mm 

Medium: 

Velocity  of  sound 

1540  m/s 

Non-lin.  parameter  (beta) 

3.5 

Attenuation 

0  dB/cm/MHz 

Number  of  harmonics 

4 

Number  of  lateral  samples 

640 

Density 

1000  kg/m3 

Maximum  depth 

100  mm 

Sampling  in  z 

0.39  mm 

Sampling  in  x  and  y 

0.308  mm 

TABLE  I 

Simulation  parameters. 


III.  Results 

The  main  simulations  are  done  for  a  set  of  sparse  array  lay¬ 
outs  that  were  developed  during  the  project  “Real-time  3D  Ul¬ 
trasound  Imaging  System  with  Advanced  Transducer  Arrays 
(NICE)”  where  the  objective  was  to  lay  the  foundation  of  a  real¬ 
time  3D  imaging  ultrasound  system.  Part  of  this  project  has 
been  done  at  this  University,  specifically  the  construction  of  the 
sparse  transducer  for  fundamental  imaging.  In  the  NICE  project 
the  geometries  were  evaluated  by  simulation  and  experiment  for 
linear  wave  propagation  [1]. 

When  using  harmonic  imaging,  the  transmitted  frequency  has 
to  be  lowered  because  the  double  frequency  has  to  be  within  the 
transducer  bandwidth.  This  has  not  been  done  in  these  simula¬ 
tions  because  the  grating  lobes  are  moved  to  other  locations  and 
it  is  not  possible  to  compare  with  conventional  imaging. 

The  specifications  used  here  are  the  same  as  those  used  in 
the  project.  The  center  frequency  is  3  MHz  and  308  /im  pitch 
(0.6  times  the  nominal  wavelength)  and  the  2D  array  has  50x50 
elements.  No  apodization  is  used,  and  the  corner  elements  are 
removed  to  approximate  a  circular  array.  The  parameters  used 
in  this  simulation  are  summarized  in  table  I.  Only  the  geometry 
is  changed  in  the  different  simulations. 
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(a)  First  harmonic. 
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(b)  Second  harmonic. 

Fig.  1.  Nonlinear  grating  lobes 

IV.  Response  from  a  point  scatterer 

A  first  order  approximation  to  calculate  the  pulse-echo  re¬ 
sponse  is  used.  It  is  based  on  the  product  of  the  continuous 
wave  of  the  transmit  and  the  receive  responses,  [6]. 

The  two-way  response  was  calculated  by  assuming  nonlin¬ 
ear  propagation  towards  a  point  scatterer,  and  linear  propagation 
back  to  the  transducer.  The  linear  propagation  of  the  reflected 
beam  is  a  valid  assumption.  The  reason  for  this  is  that  the  reflec¬ 
ted  signal  is  so  weak  that  it  generates  almost  no  distortion  of  this 
signal.  Simulation  of  the  receiving  frequency  is  set  to  the  double 
of  the  transmission  frequency  because  it  is  the  second  harmonic 
that  is  interesting  to  study.  This  way  to  calculate  the  pulse-echo 
is  valid  in  the  focus  of  the  beam. 

The  receive  and  transmit  fields  are  multiplied  and  compared 
to  a  two-way  fundamental  image. 

V.  Nonlinear  grating  lobes 

When  using  the  nonlinear  simulator,  a  new  grating  lobe  was 
found  at  the  second  harmonic.  This  new  lobe  is  a  kind  of  grating 
lobe  and  appears  at  half  of  the  sine  to  the  angle  of  the  funda¬ 
mental. 

It  is  known  that  fingers  are  created  at  the  second  harmonic. 
The  first  report  on  fingers  is  from  1973.  They  were  found  exper¬ 


imentally  by  Lockwood  et  al.  [7],  and  later  explained  theoretic¬ 
ally  by  Tj0tta  et  al.  [8],  Fingers  appear  at  the  second  harmonic 
and  they  are  the  lobes  that  appear  approximately  between  the 
sidelobe  peaks  at  the  fundamental. 

The  nonlinear  grating  lobe  phenomenon  is  illustrated  in  fig¬ 
ure  1.  The  scaling  of  the  images  is  reduced  to  -25  dB  in  the  first 
harmonic  figure  and  -50  dB  in  the  second  harmonic.  This  has 
been  done  to  make  the  grating  lobes  appear  clearer  since  noise 
is  removed.  The  lobes  in  the  middle  of  the  second  harmonic  are 
the  nonlinear  grating  lobes.  The  phenomenon  resembles  fingers. 
The  ratio  between  sine  to  the  grating  lobe  and  this  nonlinear 
grating  lobe  is  found  to  be  approximately  .49  in  this  example. 

VI.  Conclusion 

In  our  study  we  have  found  that  a  grating  lobe  at  sin  <f>g  at  the 
fundamental  fundamental  frequency  results  in  a  new  secondary 
grating  lobe  at  0.5  •  sin  <j>g  at  the  second  harmonic.  An  explana¬ 
tion  for  this  phenomenon  can  probably  be  found  from  the  theory 
of  two  intersecting  beams  in  a  nonlinear  medium. 

Although  details  have  not  been  given  here,  we  have  found 
that  sidelobes  have  been  reduced  by  a  significant  amount  for  the 
sparse  arrays  we  have  considered.  These  arrays  have  between 
208  and  880  connected  elements  on  the  receiver  and  transmitter 
out  of  the  total  2500  elements.  In  general  the  two-way  side- 
lobe  level  at  the  second  harmonic  has  been  reduced  by  20-25  dB 
compared  to  the  fundamental. 
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Abstract 

At  a  previous  symposium  [1]  we  reported  observations  of  streaming  around  a  needle  hydro¬ 
phone.  The  hydrophone  was  positioned  near  the  focus  of  a  focusing  ultrasound  transducer, 
and  streaming  was  observed  even  when  using  single  bursts.  No  streaming  was  found  without 
the  probe.  It  was  proposed  that  the  origin  of  this  streaming  could  be  exitation  of  evanescent 
waves  on  the  probe  surface.  Presently  we  report  results  from  a  new  experiment  on  this  subject. 
Streaming  was  observed  and  measured  using  tiny  tracer  particles  which  were  recorded  with  a 
CCD-video  camera  and  a  computer  frame-grabber.  The  frames  were  analysed  with  a  Matlab- 
program  for  obtaining  two-dimensional  velocity  measurements.  The  results  show  that  we  also 
observe  streaming  in  free  field,  and  even  when  using  single  bursts. 

Experimental  setup 

There  are  several  possible  ways  to  measure  acoustic  streaming  [2-4],  or  more  generally  motion 
in  fluids.  In  our  experiments  we  chose  Particle  Image  Velocimetry  (PIV)  for  several  reasons. 
PIV  experimental  setups  can  be  made  simple  and  inexpensive. 
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Figure  1:  Experimental  setup. 


The  experimental  setup  is  shown  in  Figure  1.  The  tank  has  dimensions  0.3  x  0.3  x  0.5  m, 
and  is  filled  with  tap  water  which  holds  a  temperature  of  20°  C.  The  focusing  transducer  is  a 
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single  element  of  PZT-4D,  which  has  a  nominal  resonance  frequency  of  1MHz,  a  diameter  of 
36  mm  and  a  nominal  focal  range  of  100  mm.  The  transducer  is  mounted  air  backed  in  a  brass 
housing,  and  it  is  powered  by  a  HP8116A  burst/function  generator  amplified  with  a  ENI240 
power  amplifier.  In  addition  a  HP3312A  function  generator  was  used  to  get  a  lower  burst 
repetition  frequency  than  1  Hz  (which  is  the  maximum  burst  repetition  frequency  possible 
from  the  HP8116A).  A  Saven  CI-220  laser  together  with  an  external  lens  was  used  to  produce 
a  laser  'sheet  to  illuminate  a  thin  slice  of  the  water.  The  tracer  particles  used  are  "polyamid 
seeding  particles",  PSP-50,  provided  by  Dantec  Measurement  Technology.  The  density  of  the 
particles  is  1.09  g/cm3,  and  the  mean  diameter  is  50  pm.  Unfortunately  no  data  on  sound 
velocity  in  these  particles,  or  the  material,  are  available  from  the  producers.  If  the  sound 
velocity  differs  a  great  deal  from  that  of  water,  it  is  possible  that  we  are  studying  the  direct 
effect  of  the  radiation  pressure,  and  not  the  acoustic  streaming.  However,  by  using  colour 
dye  in  the  water  one  can  make  rough  estimates  of  the  acoustic  streaming,  and  compare  this 
with  the  quantitative  measurements  when  using  the  tracer  particles.  Approximately  the  same 
velocities  are  obtained,  indicating  that  we  are  indeed  studying  acoustic  streaming.  To  capture 
the  streaming  a  miniature  black  and  white  CCD-camera  was  used,  model  MTV-261CM, 
produced  by  Elfa.  The  camera  sends  signals  to  a  frame-grabber  (ATI  All-In- Wonder  128) 
mounted  in  a  Compaq  Deskpro  2000  PC,  where  the  pictures  are  stored  for  later  processing. 
A  clock  was  used  to  synchronize  the  pulse/function  generator  and  the  frame-grabber.  As  a 
basis  for  calculating  streaming  velocities  MatPIV  version  1.4  was  used  (5).  MatPIV  contains  a 
series  of  Matlab  files  written  by  Johan  Kristian  Sveen  to  analyse  PIV  data.  Cross-correlation 
is  used  to  determine  the  displacement  in  two  captured  pictures,  two-dimensional  FFT  is  used 
to  calculate  the  correlation. 

Results 

Acoustic  streaming  was  studied  under  several  different  conditions.  The  following  results  show 
acoustic  streaming  along  the  sound-axis  of  the  transducer.  The  primary  parameters  that  were 
varied  in  this  study  are  burst  repetition  frequency,  burst  length  and  burst  amplitude.  In  the 
results  shown  in  this  paper  burst  amplitude  provided  to  the  power  amplifier  was  held  constant 
at  a  maximum  of  IV.  In  Figure  2  burst  length  was  varied  from  100  ps  to  400  [is  in  steps  of 
100  fis,  while  burst  repetition  frequency  was  held  constant  at  10  Hz. 


Figure  2:  Acoustic  streaming,  burst  repetition  frequency  10  Hz. 
a)  along  the  sound  axis,  b)  transverse  to  the  sound  axis. 

Figure  2a  shows  streaming  along  the  sound  axis,  while  Figure  2b  shows  the  transverse  profile 
of  the  streaming  at  110  mm  from  the  sound  source.  The  streaming  was  measured  after  it  had 
reached  its  maximum  value. 
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Limitations  in  the  experimental  setup  reduces  the  maximum  streaming  velocity  which 
can  be  measured.  The  main  problem  was  lack  of  computer  power,  which  limited  the  picture 
quality  and  the  ability  to  obtain  more  than  10  frames  per  second.  It  is  seen  in  Figure  2  that 
when  the  streaming  velocity  reached  10  mm/s  the  uncertainty  in  the  measurements  increased 
a  great  deal.  By  reducing  the  burst  repetition  frequency  to  1  Hz,  and  again  varying  the  burst 
length  one  obtained  the  results  shown  in  Figure  3. 


Figure  3:  Acoustic  streaming,  burst  repetition  frequency  1  Hz. 
a)  along  the  sound  axis,  b)  transverse  to  the  sound  axis. 

The  streaming  decreased  with  a  factor  5  when  the  burst  repetition  frequency  was  reduced 
from  10  Hz  to  1  Hz.  When  we  used  a  single  acoustic  burst,  and  the  length  of  the  burst  was 
varied  like  in  the  previous  cases,  one  could  still  observe  acoustic  streaming.  The  velocities 
are  shown  in  Figure  4.  The  streaming  decreased  with  a  factor  3-4  from  the  case  with  burst 
repetition  frequency  1  Hz.  But  the  streaming  was  still  observable.  In  this  velocity  range  the 
uncertainty  is  large,  but  the  results  give  an  indication  of  the  shape  and  velocity  of  the  single 
burst  streaming. 


Figure  4:  Acoustic  streaming,  single  burst, 
a)  along  the  sound  axis,  b)  transverse  to  the  sound  axis. 

The  relation  between  burst  lengths  and  streaming  velocities  for  the  1Hz  and  the  single  burst 
cases  are  shown  in  Figure  5. 
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Figure  5:  Variation  in  streaming  velocity  with  burst  length,  a)  1  Hz,  b)  single  burst. 

Decay  and  rise  times  of  the  streaming  were  also  investigated.  Figure  6a  shows  rise  time 
when  using  a  burst  repetition  frequency  of  10Hz  and  a  burst  length  of  100  (is.  The  sound 
was  turned  on  at  time  0  s,  and  turned  off  after  about  6  s.  In  Figure  6b  a  single  acoustic  burst 
with  burst  length  400  ps  was  used.  The  burst  was  transmitted  to  the  fluid  at  time  0  s.  As 
can  be  seen  in  Figure  6a,  it  took  longer  time  for  the  streaming  to  reach  its  maximum  than  in 
the  single  burst  case. 


Figure  6:  Acoustic  streaming,  decay  and  rise  time, 
a)  10  Hz  burst  repetition  frequency,  b)  single  burst. 

When  we  used  a  single  acoustic  burst  the  maximum  velocity  was  reached  right  after  the 
burst  had  passed  through  the  fluid.  In  the  10  Hz  case  if  the  sound  is  continued  after  6  s,  the 
streaming  ends  up  with  the  velocity  shown  in  Figure  2. 

Conclusion 

In  this  study  we  have  investigated  acoustic  streaming.  It  has  been  shown  that  even  a  single 
burst  can  produce  observable  acoustic  streaming.  The  streaming  reaches  its  maximum  value 
faster  in  time  with  a  single  burst  than  when  using  more  than  one  burst.  It  is  also  observed 
that  as  the  burst  repetition  frequency  is  reduced,  the  streaming  along  the  sound  axis  falls  off 
faster  after  it  has  reached  its  peak.  The  maximum  also  seems  to  move  closer  to  the  sound 
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source,  and  the  transverse  profile  broadens.  Contrary  to  our  previous  findings  this  streaming 
is  not  dependent  on  evanescent  waves,  but  exist  in  free  field. 
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ABSTRACT  Multiple -frequency  acoustic  data  collected  on  routine  surveys  may  improve  the  quality  of  scrutinised 
data  as  compared  to  single  frequency  data.  Multiple  frequency  data  is  also  used  to  investigate  herring  avoidance 
reactions  during  passage  of  a  research  vessel.  Utilising  the  strong  difference  in  directivity  of  adult  herring  at  18,  38 
and  120  kHz,  a  comparison  of  the  integrated  backscattered  echo  intensity  as  a  function  of  depth  suggests  that  most 
of  the  reduction  in  echo  energy  during  passing  must  be  due  to  density  dilution,  rather  than  to  diving. 

INTRODUCTION  Fish  is  the  second  most  important  Norwegian  merchandise,  and  its  importance  is  increasing. 
The  Institute  of  Marine  Research  (IMR)  is  responsible  for  assessment  of  important  fish  species  in  Norwegian  waters, 
and  so  far  the  use  of  acoustical  methods  as  the  first  step  in  abundance  estimation  seem  to  give  as  least  as  accurate 
estimations  for  most  species  as  any  other  method  [1],  IMR  collects  large  amounts  of  acoustic  data  during  1500  ship- 
days  at  sea  per  year,  and  these  data  needs  to  be  processed  during  the  2  -  3  hours  available  time  per  day  at  each 
cruise.  Fig.  1  illustrates  the  procedure  to  scrutinise  the  acoustic  data. 

Each  species  is  frequently  surveyed  at  the  time  of  year 
when  it  is  most  favourably  distributed.  On  acoustic 
surveys,  wily  a  few  species  are  scrutinised  to  the  highest 
accuracy  due  to  the  limited  capacity  for  biological 
sampling.  The  experience  of  the  operator  is  essential  in 
the  scrutinising  process.  Scattering  models  applied  on  the 
multiple  frequency  data  are  used  to  either  retain  the 
acoustic  returns  from  the  wanted  targets  only,  or  to 
exclude  the  returns  from  the  unwanted  targets,  and  thus, 
improve  the  quality  of  the  scrutinised  acoustic  data  as 
compared  to  data  achieved  from  single  frequency 
techniques.  Requirements  for  collection  of  optimal 
multiple  frequency  data  were  discussed  by  Komeliussen 
12]  as  well  as  some  methods  for  visualisation  of  such  data. 

Acoustic  signatures  are  used  to  discriminate  between 
simple  classes  of  acoustic  targets  [3], 

Multiple  frequency  data  can  also  be  used  to  investigate  vessel  avoidance,  which  is  regarded  the  most  important  bias 
in  acoustic  abundance  estimates  of  pelagic  fish  [4].  During  wintering,  the  important  pelagic  stock  of  the  Norwegian 
Spring  Spawning  herring  ( Clupea  Harengus  L.)  is  found  in  horizontally  concentrated,  distinct  layers,  and  favourably 
structured  fra1  mapping  and  abundance  estimation  [5],  However,  the  measured  hourly  mean  acoustic  backscattering 
coefficient  of  the  over-wintering  herring  is  3  times  larger  during  the  6  hours  of  twilight  as  compared  to  the  18  hours 
of  night.  The  reaction  by  herring  to  vessel  noise  [6]  is  according  to  modelling  theory  [7][8]  a  combination  of 
horizontal  avoidance  and  diving,  depending  on  when  the  flight  reaction  occur  in  relation  to  the  position  erf  herring 
relative  to  the  vessel.  Adult  herring,  27  -  37  cm  are  quite  directive  targets  at  the  echo  sounder  frequencies  used  [9], 
and  a  systematic  change  in  tilt  angle  distribution  may  dramatically  change  the  backscattered  energy,  even  at  constant 
within-beam  fish  density.  For  dense  layers  of  herring,  we  suggest  that  the  sharp  directivity  of  herring  may  be  used  to 
extract  information  on  the  behaviour  of  the  fish  within  the  vertical  beams  of  the  passing  vessel.  All  topics  with 
respect  to  herring  vessel  avoidance  may  not  be  investigated  from  the  measuring  platform,  but  the  relation  between 
dilution  and  diving  may  be  studied.  The  results  may  be  important  for  modelling  of  herring  avoidance  for  correction 
of  measurement  bias,  but  may  also  be  used  during  selection  of  optimal  echo  sounder  frequencies  for  herring  research. 

DATA  COLLECTION  IN  GENERAL  Calibrated  multiple  frequency  acoustic  data  [  10][  1 1]  for  both  applications 
below  were  collected  in  1999  from  R/V  “G.  O.  Sars”  cruising  at  normal  vessel  speed,  1 1  knots.  The  transducers  were 
mounted  close  together  on  a  protruding  keel  [  12].  The  transducers  at  resonance  frequencies  18,  38,  120  and  200  kHz 
have  nearly  identical  (1 1°,  7°,  7°  and  7°)  half  power  beam  angles.  A  special  software  version  of  the  Simrad  EK500 
echo  sounder  [13]  was  used  to  achieve  an  identical  burst-length  and  a  high  digital  sampling  rate  at  all  frequencies. 
The  data  was  corrected  for  noise  [14][15]  and  scrutinised  in  the  post-processing  system  [16]  after  collection. 


Fig.  1.  Interpretation  of  acoustic  data. 
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APPLICATION  I:  EXTRACTION  OF  MACKEREL  FROM  MIXED  ACOUSTIC  RECORDINGS  The  data 
was  collected  in  the  North  Sea  in  October  1999,  and  the  biological  sampling  was  optimised  for  mackerel  (Scomber 
Scombrus  L).  The  acoustic  signatures  are  integrated  into  the  post-processing  system  [16]  to  generate  a  masking 
matrix  for  extraction  of  the  acoustic  returns  from  the  target  species  mackerel  only,  and  thus,  to  exclude  echoes  from 
all  other  targets  The  masking  matrix  contains  the  number  1  for  the  elements  representing  volume-segments 
classified  as  mackerel,  and  zero  elsewhere.  Fig.  2a-d  below  shows  the  original  acoustic  data  corrected  for  noise.  The 
two  curves  in  the  upper  part  of  the  figure  shows  the  acoustic  signature  of  mackerel  (left)  and  some  other  species 
(right)  Since  we  only  had  biological  trawl-samples  for  mackerel  here,  we  are  not  able  to  tell  which  species  the  right 
signature  represents  (or  even  if  it  is  a  single  or  many  species).  Fig.  2e  shows  the  200  kHz  echogram  multiplied  by  the 
masking  matrix  and  Fig.  2d  shows  the  colour  scale  used  to  visualise  the  acoustic  data. 


s  v-f  re  q  u  en  c  y- de  p  en  d  e  nc  y  relative  to  sv3S 


Fig  2  Original  echograms  from  the  North  Sea  after  noise  removal  (a-d)  and  a  200  kHz  echogram  where  only  the 
target-category  “ mackerel ”  is  retained.  This  is  achieved  by  multiplying  the  original  200  kHz  data  with  the  masking 
matrix.  Fig.  e  shows  the  colour  scale  used. 


APPLICATION  n:  VESSEL  AVOIDANCE  FROM  MULTIPLE  FREQUENCY  DATA  The  data  was  collected 
in  Ofotfjord  January  1999  in  Northern  Norway,  in  the  wintering  area  of  the  Norwegian  Spring  Spawning  herring.  A 
typical  night  time  recording  of  herring  at  the  commonly  used  echo  sounder  frequency,  38  kHz,  is  shown  in  Fig.  3. 


Avoidance  is  expected  to  be  strong  in  the  upper 
layers  during  night.  No  significant  avoidance  is 
expected  to  occur  in  the  deeper  layers  at  night, 
below  150  m  depth,  or  in  the  daytime  recordings, 
where  all  herring  is  registered  deeper  than  150 
m.  The  deeper  layers,  both  day  and  night  will 
therefore  later  be  used  as  “reference”  layers,  with 
assumed  insignificant  vessel  avoidance.  In  these 
layers,  the  herring  is  observed  dorsally  at  three 
simultaneously  working  echo  sounder 
frequencies.  After  calibration  and  noise- 
correction  the  observed  acoustic  quantity  between 
the  depths  zi  and  z2  is  the  mean  area  acoustic 
backscattering  coefficient  with  units  [m2/n.mi.2]: 
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Fig.  3.  Two  nautical  miles  of  a  typical  night  time  recording  of  herring 
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The  effect  of  observation  volume  is  accounted  for  in  the  equivalent  beam  angle  of  the  particular  transducer  used.  The 
mean  volume  backscattering  coefficient,  <Sv>,  is  the  product  between  the  mean  volume  density  of  fish,  <fK  >,  and 
the  mean  backscattering  cross  section,  <o>,  of  the  fish:  <s ^>  =  <pvxa>  (2) 


Further,  the  mean  backscattering  cross  section,  <a>,  for  a  directive 
scatterer  is  the  convolution  of  its  directivity  pattern,  a  =  a(9),  and  its 
tilt  angle  distribution  p  =  p(0)  [15].  It  is  assumed  that  both  the 
volume  density  and  the  tilt  angle  distribution  is  the  same  in  all  three 
beams,  at  any  depth.  The  expected  differences  in  integrated  energy 
across  frequency  must  therefore  be  derived  from  differences  in  the 
fish  directivity  pattern,  o(0).  An  example  of  this  for  an  adult  herring 
is  indicated  in  Fig.  4  for  18  and  38  kllz. 

The  sharper  main  lobe  of  the  directivity  pattern  at  38  kHz  compared 
to  18  kHz,  and  even  sharper  at  120  kHz  (not  shown),  makes  the  18 
kHz  fairly  stable  for  moderate  changes  in  tilt  angle  distribution.  At 
the  two  higher  frequencies,  however,  moderate  changes,  in 
particular  in  mean  tilt  may  cause  a  severe  reduction  in  backscattered 
echo  intensity.  In  result,  if  the  fish  is  moving  from  a  natural, 
horizontal  swimming  mode  to  a  diving  mode,  it  should  affect  the 
frequencies  differently.  Actually,  by  computing  the  difference  in 
target  strength  from  the  directivity  pattern,  as  a  function  of  diving 
angle,  the  relative  drop  on  target  strength  may  be  studied  as  a 
function  of  angle.  Fig  5. 

If  the  unaffected  layers  may  be  used  as  reference  layers,  and  the  ratio 
between  integrated  echo  intensity  between  frequencies: 

ATS  =  lOlog(Ao);  Ao  =  SA.no/  Sajs  or  Ao  =  SaV  Sa,i8  (3) 

is  studied  as  a  function  of  depth,  a  substantial  drop  is  expected  in  the 
layers  where  the  fish  is  diving,  at  least  for  moderate  mean  diving 
angles.  Fig.  5.  During  data  selection,  a  detailed  scrutinising  of  one 
24-hour  cycle  of  survey  data  from  January  8  to  9  was  made  at  all 
three  frequencies,  and  stored  in  20-m  layers,  30-minute  sections. 
Further,  to  avoid  working  with  low-density  regions,  a  minimum  area 
backscattering  coefficient  of  sA  =  500  [m2/n.mi.2]  (0.0073  fish/m3) 
was  applied  as  lower  limit. 


Fig  3.  Target  strength  of  a  30 -  cm  herring 
as  a  function  of  tilt  angle  0  over  345° 


H>ing  angle  (deg) 


Fig.  5.  Expected  relative  drop  in  target 
strength  or  volume  backscattering  strength 
between  frequencies  as  a  function  of  mean 
diving  angle(smoothed  and  normalised  data). 


RESULTS  AND  DISCUSSION:  APPLICATION  1  Comparison  of  the  retained  200  kHz  data  in  Fig.  2e  with  the 
biological  trawl-samples  and  Fig.  2d  shows  a  discriminating  power  between  mackerel  and  all  other  targets.  A 
successful  generation  of  a  reasonable  masking  matrix  (which  is  extracted  from  a  target  classification  matrix)  requires 
high  quality  of  the  collected  acoustic  data  as  described  by  Komeliussen  [2].  In  reality,  it  is  not  a  single  target,  but  a 
volume-segment  that  is  classified.  Even  though  the  echo  sounder,  the  transducers  and  the  internal  transducer 
mounting  are  adapted  to  meet  these  requirements,  noise  removed,  the  data  is  not  perfect  for  the  generation  of  the 
target  classification  matrix.  Thus,  both  the  noise-corrected  multiple  frequency  acoustic  data  at  its  original  resolution 
and  a  smoothed  version  of  the  same  data  are  used  as  input  to  the  classification  system.  The  classification  system  is 
further  implemented  as  hierarchical  expert  system  where  the  first  trial  uses  a  set  of  strict  rules  for  a  volume-segment 
to  be  classified.  In  the  following  trials,  the  requirements  on  the  acoustic  data  itself  are  loosened  as  compared  to  the 
first  trial,  but  also  the  classification  of  the  neighbouring  volume-segments  are  used. 

For  mackerel  in  particular,  the  classification  seems  to  be  successful,  and  the  automatic  classification  system  also 
reduces  the  time  needed  to  scrutinise  the  acoustic  data.  For  other  targets  than  mackerel,  the  classification  seems  to  be 
promising,  but  these  can  only  partly  be  verified  by  the  biological  sampling.  This  is  due  to  the  insufficient  depth- 
resolution  and  limited  capability  of  the  present  sampling  tools  in  separating  the  different  types  of  biological  targets. 

RESULTS  AND  DISCUSSION:  APPLICATION  II  The  very  dense  concentrations  of  herring  in  the  survey  area 
may  defend  the  use  of  the  120  kHz  system  well  outside  the  normal  use  for  plankton  and  single  fish,  125  -  150  m, 
and  the  signal  to  noise  ratio  is  still  good  at  300  m.  Fig.  6  shows  the  depth  dependency  of  the  ratio  of  the  mean  area 
backscattering  coefficients,  ri=Ci(sA3g/sA.is)  and  r2=C2(sA.i2o/sA.i8)  where  the  constants,  Ci  and  c2  are  the  normalising 
constants. 
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For  diving  being  the  only  avoidance  reaction  of  herring,  ri  and  r2  should  increase  rapidly  with  depth,  at  least  down  to 
100  m  below  the  surface,  and  then  stabilise  at  a  maximum  value  at  larger  depths.  Inspection  of  Fig.  3  shows  that 
even  a  moderate  mean  diving  angle,  e.g.  5  degrees,  should  give  a  significantly  larger  effect  than  the  one  observed. 
For  density  draining  being  the  major  avoidance  reaction,  there  should  non,  or  only  a  weak  depth  dependency  of  ri 
and  r2.  The  limited  response  seen  here  therefore  suggests  that  the  avoidance  observed  in  the  shallow  layers  of  herring 
at  the  wintering  area  are  mainly  caused  by  horizontal  dilution.  The  discontinuity  at  about  100  m  is  however 
interesting,  as  ri  here  turn  back  towards  1.0  in  the  more 
shallow  layers.  This  may  indicate  either  that  the  fish  is 
swimming  more  horizontal  again,  or  that  the  fish  is  diving  at 
a  very  steep  angle,  as  suggested  in  [8].  As  seen  from  Fig.  5, 
the  ratio  also  turns  against  normal  values  at  high  diving 
angles.  The  clear  but  weak  response  seen  in  both  measures 
may  nevertheless  indicate  some  diving,  in  particular  between 
100  -  200  m  where  the  deviations  are  significant  at  both 
frequencies.  Since  the  magnitude  of  the  response  is  low,  other 
effects  may  also  affect  the  results.  These  are  depth  dependent 
target  strength,  frequency  dependent  extinction  and  others. 

Before  a  more  firm  conclusion  can  be  made,  these  factors 
should  be  investigated,  and  comparative  data  from  a  drifting 
ship,  with  stopped  main  engine  should  be  collected.  This 
would  easily  isolate  acoustic  and  behavioural  effects  on  the 
measured  ratios. 

CONCLUSIONS  Multiple  frequency  acoustic  data  may  significantly  improve  the  quality  of  the  scrutinising  process 
in  acoustic  surveys  within  the  2  -  3  hours  available  time  available  per  day.  The  multiple  frequency  data  may  also  be 
used  extract  sensible  classification  parameters,  as  well  as  investigate  the  fish  vessel  avoidance.  Thus,  some  of  the 
uncertainty  of  the  estimated  fish  stock  abundance  is  reduced. 


Fig.6.  Depth  dependency  of  the  normalised  ratio  of 
the  mean  area  acoustic  backscattering  coefficients. 
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1.  Introduction 

Multipath  ultrasonic  transit-time  meters  for  gas  and  liquid  flow  measurement  (USM)  have 
already  been  developed  to  a  stage  where  they  are  considered  as  competitive  alternatives  to  the 
more  conventional  orifice  plate,  turbine  and  positive  displacement  meters  for  fiscal  metering  of 
natural  gas  and  oil  products,  particularly  for  transmission  line  applications.  USM  technology 
offers  significant  operational  advantages  such  as  no  moving  parts,  non-intrusive  measurement  (no 
obstruction  of  flow),  no  pressure  loss,  and  bi-directional  operation  (reducing  need  for  pipework). 
Compact  metering  stations  can  be  constructed  on  basis  of  the  large  turn-down  ratio  of  USMs 
(40:1  or  larger  for  gas  meters,  20:1  or  larger  for  liquid  meters  (tentatively),  reducing  the  need  for 
a  multiplicity  of  meters  to  cover  a  wide  flow  range),  and  the  short  upstream/downstream 
requirements  with  respect  to  bends  (10D  and  5D  for  gas  meters,  typically).  Measurement 
possibilities  are  offered  which  have  not  been  available  earlier,  such  as  fast  time  response  and 
flow  monitoring  (e.g.  pulsating  flow,  flow  velocity  profile;  sound  velocity  profile),  and  self¬ 
checking  capabilities  (from  sound  velocity,  signal  level,  etc.).  In  gas  there  has  been  demonstrated 
potentials  of  additional  information  such  as  gas  density  and  calorific  value  determination.  The 
potential  of  remote  operation  of  USMs  is  an  interesting  perspective. 

For  natural  gas,  the  first  generation  of  USMs  have  been  on  the  market  for  about  5-10  years  [1-3]. 
USMs  have  demonstrated  their  capability  to  provide  metering  accuracy  within  national  regulation 
requirements  [4,5].  Better  than  ±0.7  %  uncertainty  (of  measured  value)  is  being  reported,  as 
required  for  custody  transfer  in  large  commercial  pipelines.  In  appropriate  applications,  multi- 
path  ultrasonic  meters  can  offer  significant  cost  benefits.  USM  technology  is  increasingly  gaining 
acceptance  throughout  the  industry,  and  is  today  in  use  in  gas  metering  stations  onshore  and 
offshore.  ISO  standardization  work  has  been  started  [6]. 

For  metering  of  oil  and  petroleum  products,  USMs  for  liquids  have  for  many  years  represented  a 
robust  alternative  in  non-fiscal  applications.  USMs  have  recently  been  introduced  also  for  fiscal 
metering  [7],  and  0.15  -  0.25  %  uncertainty  is  claimed,  based  on  in-situ  flow  calibration  (using 
prover).  The  petroleum  industry  is  at  present  gaining  field  experience  with  this  new  fiscal  liquid 
metering  technology  [8]. 

On  the  other  hand,  a  relatively  small  measurement  error  in  the  flow  rate  measured  by  a  USM  can 
easily  translate  into  a  large  sum  of  money.  For  example,  an  estimate  of  Norway’s  export  of 
natural  gas  in  2000  is  about  6000  million  USD.  Present-day  technology  for  fiscal  and  sales 
metering  of  gas  is  at  a  level  of  0.5-1  %  uncertainty  (of  measured  value).  A  systematic  error  in  the 
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flow  measurement  of,  say,  0.5  %  (which  is  still  within  usual  national  regulations  for  gas 
measurement  [5])  would  translate  into  a  measurement  uncertainty  corresponding  to  an  annual 
value  of  about  30  million  USD  for  this  2000  gas  export  estimate.  This  means  that  even  an 
uncertainty  of  about  0.2-0.3  %  in  some  cases  is  considered  as  uncomfortable  in  the  industry. 
Investigations  have  shown  that  there  are  un-exploited  potentials  in  USM  technology,  such  as  to 
reduce  systematic  errors  (cf.  [9]). 

For  example,  effects  of  flow  velocity  profiles  have  been  widely  studied  with  respect  to  USM 
integration  methods  (multipath  effects),  and  is  still  an  area  of  active  research  experimentally  and 
theoretically.  However,  the  effects  of  flow  velocity  profiles  on  the  acoustic  transit  times  have 
rarely  been  studied.  Such  effects  will  be  systematic,  and  can  be  of  relevance,  especially  for  gas 
meters  operating  at  high  flow  velocities.  They  will  influence  both  on  the  measured  flow  velocity 
and  sound  velocity. 

The  measured  sound  velocity  in  a  USM  has  traditionally  been  used  for  self-check  of  the  meter. 
Recently,  algorithms  have  been  developed  for  measurement  of  the  gas  density  based  on  the 
measured  sound  velocity  [10].  This  means  that  precise  measurement  of  the  sound  velocity  is  now 
of  even  higher  importance  than  a  few  years  ago. 

Even  in  relatively  “simple”  pipework  applications,  such  as  straigth  pipes  (relevant  e.g.  for 
metering  at  gas  and  liquid  transmission  networks),  the  axial  flow  velocity  profiles  may  take 
significantly  different  “shapes”,  depending  on  the  Reynolds  number  of  the  flow,  Re.  The  shape  of 
actual  profiles  influences  on  the  USM  measurement. 

In  metering  stations  where  compactness  is  important  (e.g.  offshore),  complex  installation 
conditions  (pipe  bends,  flow  conditioners,  etc.)  cause  disturbed  flow  velocity  profiles  which 
influence  on  the  USM  measurement.  This  concerns  both  axial  and  transversal  flow  velocity 
components. 

Transversal  flow  components  will  occur  especially  when  the  USM  is  installed  downstream  bends 
and  other  obstructions  of  the  pipe  flow.  It  is  expected  that  downstream  a  double  bend  out  of 
plane,  the  transversal  flow  regime  is  typically  a  swirl,  while  downstream  a  single  bend,  a  cross 
flow  is  typically  established.  Such  transversal  flow  components  can  contribute  to  the  flow 
measurement  performed  by  each  acoustic  path.  There  are  examples  from  installations 
downstream  double  bends  out  of  plane  where  the  transversal  flow  components  can  be  10  %  of  the 
axial  flow  component,  or  larger.  In  practice,  the  transverse  flow  components  will  often  neither  be 
a  symmetric  swirl  nor  a  symmetric  cross  flow,  but  instead  some  kind  of  an  asymmetric  variant  of 
either  swirl  or  cross  flow,  or  something  inbetween. 

In  the  present  work,  the  accuracy  of  the  traditional  expressions  used  in  today’s  USMs  for 
measurement  of  flow  velocity  and  sound  velocity,  with  respect  to  flow  profile  effects  on  transit 
times,  is  investigated  by  use  of  ray  theory.  The  influences  of  axial  and  transversal  flow  profiles 
are  both  studied.  Computational  fluid  dynamics  (CFD)  modelling  of  pipe  flow  is  used  to 
establish  representative  analytical  expressions  for  acoustic  path  profiles  (axial  and  transversal). 
The  analytical  path  profiles  are  then  used  as  input  to  a  ray  theory  model  for  sound  propagation  in 
pipe  flow,  to  evaluate  effects  of  more  realistic,  non-uniform  flow  profiles  (axial  and  transversal) 
on  the  acoustic  transit  times.  The  ray  theory  simulation  results  are  compared  with  the  simpler 
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and  more  traditional  analytical  expressions  for  the  acoustic  transit  times,  derived  for  the 
simplified  case  of  uniform  axial  and  transversal  flow  profiles.  The  resulting  errors  in  flow 
velocity  and  sound  velocity  for  the  traditional  USM  functional  relationships  are  evaluated  and 
discussed. 


2.  Current  theoretical  basis  for  USM  flow  metering  methodology 

A  USM  measures  the  axial  (x)  component  of  the  volumetric  flow  rate  at  line  conditions  (with 
respect  to  pressure,  temperature  and  fluid  quality),  qv,  defined  as 


\\vx{y,z)dydz  , 


A 


(1) 


where  A  is  the  cross  sectional  area  (in  they;,  z  -  plane)  of  the  pipe,  cf.  Fig.  1,  and  v,  is  the  axial  (x) 
component  of  the  flow  velocity.  For  circular  cross-section,  the  double  integral  in  Eq.  (1)  can  be 
written  as  a  single  integral  [11,9], 


R  -jR2-y2  r 

=  J  J  v*  O'. z)dzdy  =  2  J  Jit2  -y2vx  (y)dy ,  (2) 

-R  -R 


where  R  =  D/2  is  the  inner  radius  of  the  pipe  and 


vxOO  = 


(3) 


is  the  average  axial  flow  velocity  (the  line  integral)  over  the  chord  with  lateral  position  y. 


Fig.  1.  Schematic  illustration  of  a  single  path  (no.  i)  in  a  multipath  ultrasonic  gas  flow  meter  (for  downstream 
sound  propagation).  (Left:  centre  path  example  (y,  =  0);  Right:  path  at  lateral  chord  position  y,.) 


In  USMs  used  for  fiscal  metering,  the  flow  velocity  is  measured  using  a  finite  number  of  acoustic 
paths,  typically  4-6,  each  path  corresponding  to  a  chord  over  the  pipe  cross  section,  cf.  Fig.  1  [1 1, 
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9].  For  each  path,  transit  times  are  measured  electronically  for  high-frequency  ultrasonic  pulses 
propagating  across  the  pipe,  at  an  angle,  <t>n  with  respect  to  the  pipe  axis,  downstream  with  the 
flow,  and  upstream  against  the  flow.  Ultrasonic  transducers  are  used  to  transmit  and  receive  the 
signals,  and  the  measured  transit  times  are  corrected  for  time  delays  due  to  electronics,  cables, 
transducers,  diffraction,  transducer  cavities,  etc.  For  each  acoustic  path,  the  difference  between 
the  upstream  and  the  downstream  propagating  transit  times  is  proportional  to  the  average  flow 
velocity  along  the  acoustic  path.  Multiple  acoustic  paths  are  used  to  sample  the  flow  velocity 
profile  in  the  pipe  at  a  set  of  discrete  chords,  to  improve  the  metering  accuracy. 

In  this  approach,  thus,  Eq.  (2)  is  solved  by  numerical  integration,  giving 


(4) 


where  Nc  is  the  number  of  chords  (typically  4-5),  w)  are  the  integration  weight  factors,  j  =  1,  ..., 
Nc,  V;  r  =  v_ (v )  and  chord  no.  j  is  located  at y  =  yj. 

Jy*  *  J 


In  practice,  the  flow  velocity  is  measured  over  N  acoustic  paths,  having  inclination  angles  to  the 
axial  flow  direction  of  typically  40°  to  50°,  and  where  the  projection  of  one  such  acoustic  path 
constitutes  one  of  the  chords  used  above.  For  each  such  acoustic  path,  the  upstream  and 
downstream  transit  times  are  measured. 


2.1  Traditional  approach;  uniform  axial  flow  and  no  transversal  flow 

For  the  simplified  case  where  the  flow  velocity  profile  is  assumed  to  be  uniform  and  purely  axial 
(i.e.  no  transversal  flow  components),  the  transit  times  of  path  no.  i  can  be  found  by  a  simple 
geometrical  approach  as  (see  Fig.  2a) 


where  L,  is  the  interrogation  length,  the  inclination  angle  of  the  acoustic  path,  and  tu  and  t2t  are 
the  upstream  and  downstream  transit  times,  respectively.  vix  is  the  average  axial  flow  velocity 

over  acoustic  path  no  /,  and  c  is  the  velocity  of  sound  (VOS).  From  Eqs.  (5)  it  is  easily  seen  that 

v.  =  hklLlhil'  (6) 

2tut2i  cos  (fii 


Alternatively,  a  more  sophisticated  and  accurate  ray  tracing  approach  can  be  taken,  as  reported  by 
McCartney  et  al  [12],  leading  to 


Eq.  (7)  is  claimed  by  McCartney  et  al.  to  be  valid  also  for  non-uniform  flow  velocity  profiles. 
However,  an  underlying  assumption  in  their  analysis  is  that  the  rays  are  straight  lines.  This  is 
possible  for  uniform  flow  profiles  only,  and  Eq.  (7)  is  therefore  valid  only  for  uniform  axial  flow 
profile  and  no  transversal  flow  components  [11,9]. 

It  is  interesting  to  note  that  Eq.  (7)  leads  to  exactly  the  same  expression,  given  by  Eq.  (6),  for  the 
average  flow  velocity  over  the  acoustic  path,  vix,  as  the  simplified  geometrical  approach 
described  above  and  illustrated  in  Fig.  2a. 


Fig.  2.  Schematic  illustration  of  two  simplified  “geometrical”  approaches  which  may  be  used  to  account  for  flow 
velocity  components  in  the  USM  functional  relationships: 

(a)  Traditional  approach;  uniform  axial  flow  profile  only,  and  no  transversal  flow, 

(b)  Slightly  improved  approach;  uniform  axial  flow  profile  and  uniform  transversal  flow  profile. 

Eqs.  (7)  also  leads  to  the  well-known  expression  for  the  sound  velocity, 

c  =  Wfry +  hi  Y  cos2  <t>i  +  (hi  -  hi  Y  sin2  fa 

XnhiCOsfr  •  (8) 

The  velocity  of  sound  has  traditionally  been  used  for  self-checking  in  USMs.  In  addition,  recent 
developments  have  shown  that  it  can  be  used  as  a  basis  for  calculation  of  density  or  calorific 
value  of  a  natural  gas  [10].  Through  the  velocity  of  sound,  therefore,  a  USM  can  be  used  as  a 
mass  flow  meter  or  an  energy  flow  meter,  in  addition  to  its  more  traditional  use  as  a  volumetric 
flow  meter.  In  such  new  applications,  a  sound  velocity  measurement  with  a  uncertainty  of  about 
±0.25  m/s  or  better,  will  be  needed  for  fiscal  measurement  of  density  or  calorific  value  [19].  In 
non-fiscal  applications,  a  less  accurate  sound  velocity  measurement  may  be  acceptable. 

2.2  Improved  approach;  uniform  axial  and  transversal  flow  profiles 

In  a  real  flow  metering  situation,  there  will  be  transversal  flow  velocity  components  in  addition  to 
the  axial  flow  velocity  component.  Such  transversal  flow  velocity  components  may  influence  on 
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the  transit  times,  and  thus  also  on  the  estimated  flow  and  sound  velocities.  In  the  simplest 
approximation,  both  the  axial  and  the  transversal  flow  velocity  components  are  considered  to  be 
uniform.  A  simple,  geometrical  approach  can  again  be  taken,  giving  the  following  upstream  and 
downstream  transit  times  (see  Fig.  2b): 


tu  = 


A  =- 


L 


u  c-  vix  cos  (f)i  -  vi  z  sin  fa  ’  2‘  c  +  vix  cos  fa  +  viz  sin  fa  ’ 


(9) 


where  v;  z  is  the  average  transversal  flow  velocity  component  in  the  z  direction.  From  Eqs.  (9),  c 
can  be  eliminated,  giving 


v,'X+tan^v,.z  = 


Rj  fyn  *2  i ) 
2 tut2i  cos  fa 


(10) 


as  an  improvement  relative  to  Eq.  (6).  In  actual  fiscal  flow  meters,  Eqs.  (6)  or  (10)  (for  flow 
velocity)  and  Eq.  (8)  (for  sound  velocity)  are  the  expressions  used  to  obtain  the  average  flow 
velocity  and  sound  velocity  at  acoustic  path  no.  i. 

Note  that  Eq.  (8)  is  “exact”  (within  the  ray  approximation)  for  the  uniform  axial  flow  profiles 
when  the  transversal  flow  component,  vz,  is  zero.  Eq.  (10)  is  expected  to  be  a  good 
approximation  for  uniform  axial  and  transversal  flow  profiles. 

In  Eq.  (4),  the  average  axial  flow  velocity  over  the  chord,  vcj  x ,  is  needed,  while  in  Eqs.  (6)  and 
(10)  the  average  axial  flow  velocity  over  the  inclined  acoustic  path,  vi  x,  is  involved.  It  is  usual  to 
assume  that  these  two  quantities  are  approximately  equal, 


v .  « V-  , 

JtX  l, X  * 


(11) 


for  corresponding  chord  and  acoustic  paths.  Hence,  in  this  approximation  it  is  assumed  that  the 
axial  flow  velocity  profile  is  constant  over  the  metering  volume  of  the  USM.  However,  in  order 
to  obtain  vix,  a  complication  is  that  in  the  USM,  it  is  found  in  a  linear  combination  with  the 

transversal  flow  component,  vi  z,  see  Eq.  (10),  and  not  isolated  as  needed  in  Eq.  (4).  In  order  to 
treat  this  complication,  Eq.  (4)  is  written  as 


qv  WjVj,x  *  xR2 X v/  ’ 

y=i  m 


(12) 


where  the  inclination  angles  $  and  the  integration  weights  vv;  are  chosen  to  reduce  the  influence 
of  the  transversal  flow  components.  N  is  here  the  number  of  acoustic  paths  in  the  flow  meter, 
typically  4-6,  and 
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(13) 


V.  =  * 2i  ) 

2t nt 2i  cos  </>t  ’ 
cf.  Eqs.  (6)  and  (10). 

In  the  case  where  there  is  one  acoustic  path  per  chord,  so  that  N  =  Nc,  the  weights  w,  are  equal  to 
the  weights  wf.  In  the  case  where  there  are  2  acoustic  paths  per  chords  (for  all  or  some  of  the 
chords),  the  relation  between  the  weights  wh  /  =  1,  ....  N  and  wf,  j  =  1,  ...,  Nc  will  be  more 
complex.  It  should  also  be  noted  that  it  is  only  when  there  are  2  acoustic  paths  per  chord  for 
every  chord,  i.e.  N=  2NC,  that  the  transversal  flow  velocity  component  viz  in  Eq.  (10)  will  be 

cancelled  for  any  transversal  flow  profile.  For  a  smaller  number  of  acoustic  paths,  only  certain 
transversal  flow  velocity  profiles  will  in  general  be  cancelled.  However,  in  many  cases,  it  is 
possible  to  design  the  integration  weights  so  that  the  transversal  flow  velocity  components  of 
main  interest  in  ultrasonic  flow  metering  will  be  reduced. 

The  traditional  expressions  on  which  today’s  USMs  are  based,  given  by  Eqs.  (6)  or  (10)  (for  flow 
velocity)  and  Eq.  (8)  (for  sound  velocity),  are  based  on  assumptions  of  uniform  axial  and 
transversal  flow  velocity  profiles,  as  explained  above.  In  a  real  flow  metering  situation,  the  axial 
and  transversal  flow  velocity  components  will  rarely  be  uniform,  cf.  Section  1,  and  the  traditional 
expressions  represent  approximations.  In  the  following,  the  accuracies  of  the  traditional 
expressions  Eq.  (8)  and  (10)  are  investigated,  when  more  realistic  disturbed  (non-uniform)  axial 
and  transversal  flow  profiles  are  taken  into  account. 


3.  Flow  velocity  profiles  (axial  and  transversal) 

In  Section  5,  the  traditional  USM  expressions  given  by  Eqs.  (8)  and  (10)  will  be  investigated 
using  a  ray  theory  accounting  for  non-uniform  axial  and  transversal  flow  velocity  profiles.  As 
input  to  the  numerical  ray-tracing  model,  3-dimensional  flow  profile  data  are  needed.  This  could 
in  principle  be  treated  using  CFD-calculated  flow  profile  data  directly  into  the  ray-tracing  model. 
In  the  present  work,  a  slightly  more  simplified  approach  is  used,  in  which  the  axial  and 
transversal  flow  profiles  are  described  by  analytical  formulas.  In  order  to  make  representative 
analytical  input  flow  velocity  profiles,  CFD  calculated  flow  profiles  have  been  studied  using 
CMR’s  CFD-code  MUSIC  [13,  14],  for  various  pipework  /  bend  configurations.  Depending  on 
the  installation  conditions  of  the  USM,  the  flow  velocity  profile  will  change.  The  following 
installation  conditions  have  been  used  in  the  present  study  of  CFD-calculated  flow  profiles: 

•  USM  installed  in  long  straight  pipe,  for  various  Reynolds  numbers,  see  Fig.  3, 

•  USM  installed  10D  downstream  a  single  90°  bend,  see  Fig.  4, 

•  USM  installed  10D  downstream  a  double  90°  bend  out  of  plane,  see  Fig.  4. 

The  calculated  axial  and  transversal  flow  velocity  profiles  have  been  processed  further  using 
CMR’s  uncertainty  program  for  USMs,  GARUSO,  by  use  of  3-dimensional  cubic  spline 
interpolation  of  the  3-dimensional  flow  velocity  data.  In  this  way,  the  axial  and  transversal  flow 
velocity  profiles  over  the  acoustic  paths  are  calculated  from  the  profiles  over  the  pipe  cross 
section.  As  an  example,  the  flow  velocity  profiles  (axial  and  transversal)  have  been  calculated  for 
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an  acoustic  path  configuration  which  is  relevant  for  a  commercial  USM  for  gas.  However,  the 
specific  path  configuration  used  is  not  very  important  for  the  results  of  the  study. 
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Fig.  3.  CFD-calculated  flow  velocity  profiles  for  flow  through  a  straight  pipe  at  various  Reynolds  numbers.  Based 
on  these  CFD  data,  the  flow  velocity  profile  over  a  centre  path  in  a  USM  has  been  calculated  by  a  3- 
dimensional  cubic  spline  interpolation,  using  the  GARUSO  USM  uncertainty  model  and  simulator  [11,  9, 
15], 
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Numerical  interpolation  (GARUSO) 


Fig.  4.  CFD-calculated  flow  velocity  profiles  downstream  a  single  90°  bend,  and  a  double  90°  bend  out  of  plane. 
Based  on  these  CFD  data,  the  flow  velocity  profiles  over  four  acoustic  paths  of  relevance  in  a  USM  are 
calculated  by  a  3-dimensional  cubic  spline  interpolation,  using  the  GARUSO  USM  uncertainty  model  and 
simulator  [11,  9,  15].  For  each  bend  configuration,  two  orientations  of  the  USM  relative  to  the  pipework  are 
shown  (“horizontal”  and  “vertical”  chords).  For  each  orientation  and  bend  type,  axial  (left)  and  transversal 
(right)  path  profiles  are  shown. 


Numerical  interpolation  (GARUSO) 


For  flow  in  a  straight  pipe,  Fig.  3  shows  the  3-dimensional  CFD  calculated  profiles,  and  the  flow 
velocity  profile  over  an  acoustic  path  through  the  centre  of  the  pipe.  In  such  flow,  the  transversal 
flow  velocity  will  vanish,  and  the  axial  flow  velocity  is  symmetric.  The  well-known  dependency 
of  the  profile  flatness  on  the  Reynolds  number,  which  is  demonstrated  in  the  figure,  means  that  a 
variety  of  different  flow  profiles  have  to  be  taken  into  account  in  the  discussion  of  flow  profile 
effects  on  transit  times,  flow  velocity  and  sound  velocity,  even  for  the  "simple"  case  of  a  straight 
pipe. 

For  USM  installation  10D  downstream  a  single  90°  bend,  and  a  double  90°  bend  out  of  plane, 
Fig.  4  shows  the  3 -dimensional  CFD  calculated  profiles,  and  the  flow  velocity  profile  over  four 
acoustic  paths  of  relevance  in  the  USM  referred  to  above.  Path  profiles  for  two  orientations  of 
the  USM  relative  to  the  pipework  are  shown.  In  this  case,  both  axial  (left)  and  transversal  (right) 
flow  profiles  over  the  paths  have  been  calculated.  It  can  be  seen  that  the  axial  flow  velocity 
profiles  over  the  acoustic  paths  are  asymmetric  at  some  of  the  paths.  The  spread  of  transversal 
flow  velocity  component  over  the  various  paths  is  larger  than  for  the  axial  flow  velocity  profiles. 
Here,  a  sinusoidal  profile  is  seen  in  some  cases,  and  in  other  cases,  the  profile  is  more  similar  to 
the  axial  flow  profiles. 

In  Section  5,  these  results  are  used  to  derive  typical  and  reasonable  representative  analytical  path 
profiles  (axial  and  transversal)  which  are  used  as  input  to  the  numerical  ray  model  discussed  in 
Section  4. 

4.  Ray  theory  model  for  sound  propagation  in  pipe  flow 

A  ray  theory  model  has  been  derived,  enabling  calculation  of  the  transit  time  along  acoustic  path 
no. /  =  1,  N.  In  this  theory,  non-uniform  axial  and  transversal  flow  velocity  components  can 
be  accounted  for,  and  thus  the  effect  of  non-uniform  flow  profiles  on  the  transit  times  can  be 
studied.  The  model  is  a  further  development  of  earlier  work  by  the  authors,  cf.  [16,  1 1]. 

The  model  is  based  on  the  general  ray  tracing  equations  as  given  by  Pierce  [17]: 

ds  _  1-v-s  . 

— =■  =  —s  ■  Vv  —  5  x  V  x  v - = — Vc,  (14) 

dt  c 


dx  c2  s 

—  =  v  + - — , 

dt  1-v-s 


(15) 


where  s  is  the  wave  slowness  vector,  t  is  the  time,  x(t)  describes  the  ray  trajectory,  and  v(x)  = 
(vx(x,y,z),vy(x,y,z),vz(x,y,z))  is  the  flow  velocity  vector.  In  addition  to  the  ray  approximation, 
which  is  basically  a  high-frequency  approximation  where  diffraction  effects  are  neglected,  the 
following  assumptions  are  made  (cf.  Fig.  5): 

•  Constant  velocity  of  sound,  c, 

•  2-dimensional  flow:  vy  =  0, 

•  vx  and  vz  are  independent  of  x  and  y:  vx  =  vx(z),  vz  =  vz(z). 
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Rx 


Tx 


Fig.  5.  Principle  sketch  of  the  ray  path  from  the  transmitting  transducer  (Tx)  to  the  receiving  transducer  (Rx),  for 
downstream  propagation  in  acoustic  path  no.  i. 

Eqs.  (14)  and  (15)  then  reduce  to  the  set  of  equations 


dsz  _  dvx  dvz 
dt  $x  dz  $z  dz 

~  =  vx(z)  + - - 

dt  -vz(z)s2 

dz  c2sr 

—  =  V2(z)  + - 5 - 

dt  l~vAz>x-vAz)sz 


(16) 


In  order  to  solve  Eq.  (16),  initial  conditions  are  specified  at  t  =  0.  For  downstream  propagation, 
these  initial  conditions  are,  for  x  =  0  and  z  =  -R, 


_ _ cos  ^9 _ 

c  +  vx  (-R)  cos  (p  +  vz  (-R)  sin  (p  ’ 
sin^? 

sz  - - - - . 

c  +  vx  (-R)  cos  <p  +  vz  (-R)  sin  <p 

Eq.  (16),  with  initial  conditions  given  by  Eq.  (17),  has  been  solved  using  a  fourth  order  Runge 
Kutta  method.  Analytical  expressions  for  vx  and  vz  (axial  and  transversal  flow  profiles, 
respectively)  have  been  derived  and  used,  cf.  Section  5.  An  iteration  procedure  for  <p  based  on 
Newton's  method  has  been  used  to  ensure  that  the  ray  ends  at  the  specified  receiver  point. 

The  simulation  of  the  upstream  sound  propagation  is  similar  to  the  downstream  propagation,  with 
modified  initial  conditions  (expressions  not  given  here). 

It  may  be  noted  that  Eqs.  (16)-(17)  could  have  been  obtained  also  if  Lighthill’s  [18]  set  of  basic 
ray  theory  equations  were  used  instead  of  Eqs.(14)-(15).  Note  that  the  present  solution  gives  only 
ray  transit  times,  not  the  ray  amplitude. 
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5.  Results  using  the  ray  theory  model 

Based  on  the  CFD-calculated  axial  and  transversal  flow  profiles  for  straight  pipes  and  pipes  with 
bends,  discussed  in  Section  3,  typical  but  simplified  analytical  axial  and  transversal  flow  profiles 
have  been  selected  for  use  in  the  numerical  ray  simulations.  These  flow  profiles  are  analytical 
simplifications  of  the  CFD  calculated  flow  profiles  discussed  in  Section  3.  The  following 
analytical  flow  profiles  have  been  used: 

Axial  flow  profiles: 

•  Symmetric  parabolic  flow  profile  (Laminar  flow,  Re  <  1000). 

•  Symmetric  turbulent  flow  profile  (Turbulent  flow,  Re  >  3000). 

•  Asymmetric,  turbulent  flow  profile. 

•  Uniform  flow  profile,  Re  =  oo. 

Transversal  flow  profiles: 

•  0  %  uniform  profile  (=  0). 

•  10  %  uniform  profile,  selected  for  comparison  between  simple  and  more  “real”  flow 
velocity  profiles.  (The  transversal  flow  velocity  is  10  %  of  the  average  axial  flow 
velocity). 

•  10  %  symmetric  flow  profile  (of  relevance  in  cross  flow  and  swirl  situations). 

•  Sinusoidal  flow  profile  (of  relevance  in  cross  flow  situations). 

The  selected  axial  and  transversal  flow  velocity  profiles  are  shown  in  Fig.  6. 


Fig.  6.  The  selected  axial  (left)  and  transversal  (right)  flow  velocity  profiles  used  in  the  ray  theory  simulations.  The 
axial  flow  velocity  profiles  are  scaled  to  have  average  flow  velocity  of  1,  and  the  transversal  flow  velocity  is 
scaled  to  the  average  axial  flow  velocity. 

The  following  acoustic  path  parameters  have  been  used  in  the  simulations: 

•  Pipe  diameter:  D  =  0.2  m  (»  8"). 

•  Velocity  of  sound:  c  -  400  m/s. 

•  Inclination  angle:  ^=45°. 
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Fig.  7  Calculated  acoustic  ray  paths  for  upstream  and  downstream  propagation,  for  two  types  of  axial  flow 
velocity  profiles.  Here,  <f>  =  45°,  c  =  400  m/s,  D  =  20  cm,  and  vx  =  50  m/s.  The  deviation  of  the  ray  paths 
from  a  straight  line  is  caused  by  the  non-uniformity  of  the  axial  flow  velocity  profile.  Left:  T-aminar 
(parabolic)  axial  flow  profile  (Re  <  1000).  Right:  Turbulent  axial  flow  profile  (Re  >  3000,  typically). 
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Fig.  8  The  effect  of  non-uniform  flow  velocity  profiles  on  the  estimate  of  the  average  flow  velocity  over  the 
acoustic  path  (left)  and  the  velocity  of  sound  (right).  In  the  two  upper  figures,  the  transversal  flow  velocity 
component  is  zero  and  the  effect  of  various  axial  flow  velocity  profiles  is  studied.  In  the  two  bottom  figures, 
the  axial  flow  velocity  profile  is  uniform  and  the  effects  of  various  transversal  flow  velocity  profiles  are 
studied. 


If  the  axial  flow  velocity  profile  is  uniform,  the  ray  path  between  the  transmitter  and  the  receiver 
will  be  the  straight  line.  Non-uniformity  in  the  axial  flow  velocity  profile  will  cause  the  ray  path 
to  deviate  from  this  straight  line.  In  Fig.  7,  the  ray  paths  are  shown  for  a  flow  velocity  of  50  m/s, 
for  a  laminar  (parabolic)  axial  flow  profile  (left)  and  a  turbulent  flow  profile  (right).  Transversal 
flow  velocities  are  not  taken  into  account  in  Fig.  7.  In  the  case  of  the  laminar  flow  velocity 
profile,  the  deviation  of  the  ray  path  from  the  straight  line  is  larger  than  for  the  turbulent  flow 
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velocity  profile.  This  is  because  the  turbulent  flow  velocity  profile  is  closer  to  a  uniform  flow 
profile  (which  gives  a  straight  line)  than  the  laminar  profile. 

The  deviation  of  the  ray  path  from  a  straight  line  will  also  influence  on  the  upstream  and 
downstream  transit  times  of  the  path.  This  will  again  affect  the  flow  velocity  and  sound  velocity 
measured  over  the  acoustic  path,  as  shown  in  Fig.  8.  The  figure  shows  the  deviation  from  Eq. 
(10)  (left  part)  and  Eq.  (8)  (right  part)  (which  are  used  in  flow  meters  today),  when  the  transit 
times  are  calculated  by  ray  theory  using  the  non-uniform  flow  velocity  profiles  as  input. 

In  the  upper  left  part  of  Fig.  8,  the  deviation  is  shown  for  three  different  and  purely  axial  flow 
profiles  (no  transversal  flow  present),  for  flow  velocities  up  to  50  m/s.  Today,  multipath  USMs 
for  gas  typically  operate  up  to  velocities  in  the  range  30-40  m/s.  This  upper  limit  is  continuously 
being  pushed  upwards,  driven  by  market  needs.  It  can  be  seen  in  the  figure  that  the  turbulent  flow 
velocity  profiles  studied  here  may  give  a  systematic  deviation  of  about  0.1  %  for  flow  velocities 
around  40  m/s.  This  means  that  the  effect  is  of  relevance  when  the  total  uncertainty  of  the  meter 
should  be  around  0.5-0.7  %.  The  parabolic  profile  gives  larger  deviation  to  the  uniform  flow 
profile  solution.  However,  the  parabolic  flow  velocity  profile  is  of  most  relevance  at  liquid  flow 
metering  of  low  flow  velocities.  Therefore,  this  deviation  curve  has  been  dotted  for  high  flow 
velocities. 

In  the  lower  left  part  of  Fig.  8,  the  effects  of  transversal  flow  velocity  profiles  are  shown.  In  this 
case,  the  axial  flow  velocity  profile  is  uniform,  so  that  the  effect  of  the  transversal  flow  profiles 
can  be  studied  separately.  The  flow  velocity  calculated  from  the  transit  times  is  compared  to 
v.  =  v. .x  +  tan  (j)1vl  z  which  is  obtained  for  uniform  axial  and  transversal  flow  velocity  profiles,  cf. 

Eq.  (10).  It  can  then  be  seen  that  there  are  no  significant  effect  of  the  non-uniformity  of  the 
transversal  flow  profile,  for  the  selected  flow  profiles.  This  means  that  formulas  assuming 
uniform  transversal  flow  velocity  profiles  seem  to  be  applicable  also  for  typical  flow  velocity 
profiles  in  ultrasonic  flow  metering. 

On  the  two  right  parts  of  Fig.  8,  the  flow  profile  effect  on  the  estimated  velocity  of  sound  is 
calculated.  In  the  upper  right  figure,  transversal  flow  velocity  components  are  set  to  zero,  and 
therefore,  the  effect  of  varying  axial  flow  velocity  profiles  are  studied  isolated.  It  can  be  seen  that 
the  selected  turbulent  flow  velocity  profile  can  give  an  error  in  the  velocity  of  sound  of  up  to  0.1 
m/s  for  high  flow  velocities.  An  error  of  this  magnitude  is  significant  when  applying  the  velocity 
of  sound  as  input  to  density  or  calorific  value  estimation  algorithms  [10]. 

The  lower  right  part  of  Fig.  8  shows  the  influence  of  transversal  flow  velocity  components  on  the 
sound  velocity.  The  axial  flow  velocity  component  is  here  set  to  uniform  in  order  to  be  able  to 
study  the  effect  of  transversal  flow  velocity  components  isolated.  It  can  be  seen  that  the 
appearance  of  a  10  %  transversal  flow  velocity  component  (relative  to  the  axial  flow  velocity 
component)  will  give  an  error  in  the  estimate  of  the  velocity  of  sound  of  around  0.5  m/s  or  more 
for  large  flow  velocities.  This  effect  can  be  described  by  an  analytical  formula  (not  shown  here). 
There  is  not  observed  significant  differences  in  the  error  of  the  sound  velocity  estimate  from  non- 
uniform  transversal  flow  velocity  profiles,  relative  to  the  error  of  the  sound  velocity  estimate  due 
to  the  uniform  transversal  flow  velocity  profile,  for  the  same  average  transversal  flow  velocity 
over  the  acoustic  path. 


95 


6.  Conclusions  and  future  perspectives 


A  ray  theory  model  has  been  developed  and  used  to  study  the  accuracy  of  the  traditional 
expressions  in  current  use  for  measurement  of  the  average  flow  velocity  and  sound  velocity  at 
each  acoustic  path  in  an  ultrasonic  flow  meter  (USM).  The  ray  tracing  method  has  in  addition  the 
potential  of  improving  these  methods  in  current  use. 

With  respect  to  axial  flow  velocity  profiles,  the  effects  of  non-uniformity  in  such  profiles  are 
shown  to  be  relatively  small  for  low  flow  velocities.  In  such  applications,  the  traditional  USM 
expressions  for  flow  velocity  and  sound  velocity,  Eqs.  (10)  and  (8),  respectively,  should  be 
sufficient.  However,  at  higher  flow  velocities,  typically  30  m/s  and  higher  (Mach  number  > 
0.08),  the  results  indicate  that  the  axial  flow  profile  effects  may  be  more  significant,  giving  a 
systematic  effect. 

With  respect  to  transversal  flow  velocity  profiles,  no  significant  effects  of  non-uniformity  of  such 
profiles  have  been  found.  Transversal  flow  velocity  profiles  may  then  be  treated  as  being 
uniform.  For  this  simplified  case,  improved  expressions  for  the  flow  velocity  and  the  velocity  of 
sound  (relative  to  the  traditional  expressions)  have  been  derived  (not  given  here). 

Further  work  in  this  field  should  include  a  more  extensive  investigation  with  respect  to  flow 
velocity  profiles.  For  instance,  CFD-calculated  flow  profiles  should  be  used  directly  into  the  ray 
theory  model,  as  an  improvement  relative  to  the  analytical  profiles  used  here.  The  acoustic 
propagation  model  should  preferably  be  extended  to  a  3-dimensional  description  of  flow  velocity 
(and  not  2-dimensional  as  in  the  present  work).  Probably  more  important:  -  the  flow  profile 
effects  in  straight  pipes  should  be  analysed  as  a  function  of  Reynolds  number  (and  not  solely  as  a 
function  of  flow  velocity  as  made  here,  cf.  Fig.  8). 

The  ray  tracing  approach  used  here  is  restricted,  due  to  the  high-frequency  approximation 
inherent  in  that  methodology  (cf.  Section  4),  giving  an  infinitely  thin  acoustic  ray,  and  no 
description  of  diffraction  and  beam  effects.  To  fully  evaluate  the  systematic  transit  time  effects 
discussed  here,  and  other  systematic  effects  (cf.  e.g.  [9,  15]),  a  more  comprehensive  analysis 
based  on  wave  theory  will  be  needed,  accounting  for  acoustic  diffraction  effects,  finite  beam 
interaction  with  the  flow,  etc. 
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There  arc  many  instances  in  acoustics,  scismics  and  electromagnetics  when  it  is  useful  to  represent  wave 
fields,  including  diffraction,  in  terms  of  traveling-wave  or  directional  expansions.  Field  components 
with  a  well  defined  angle-time  arrival  structure  are  produced  in  this  way.  For  definiteness,  we  consider 
mono-frequency  acoustic  wave  propagation  with  time  dependence  e-*"*,  u  being  angular  frequency. 
Denoting  horizontal  slowness  by  p, 

F(r,  z)  =  ^  J  G(z;p)H^1)(o>pr)w2p</p  (1) 

synthetizes  the  pressure  field  at  range  v  and  depth  z  in  a  laterally  homogeneous  medium  from  its 
wavenumber  components  [1],  By  introducing  a  specification  depth  where  the  medium  is  artificially 
considered  to  be  locally  homogeneous,  cf.  [2,  eh.  6]  and  [3,  sec.  1],  we  may  introduce  p-dependent 
reflection  coefficients  74  and  7 b  for  upward  and  downward  going  plane  waves,  respectively.  The 
Green’s  function  G(z)  in  (1)  can  be  written 

G(z)  =  (1-  ran)"1  S>,0(z)  =  Gi,j(z)  (2) 

;=i  »=i  j= 0 

where  the  Gij(z)  =  (7^7b)jG;iU(z)  are  component  Green’s  functions. 

A  standing-wave  (normal  mode)  expansion  appears  as  the  residue  contributions  from  the  p  such  that 
7A1B  =  1  (e.g.,  [3]),  whereas  a  traveling- wave  expansion  [4]  appears  by  applying  the  geometric  scries  for 
(1- 7.47b)-1  as  indicated.  The  index  i  is  intended  to  separate  upgoing  and  downgoing  waves  at  source 
as  well  as  receiver,  whereas  the  index  j  is  thought  of  as  a  cycle  number.  Tire  implied  decomposition 
of  P(r,  z)  becomes 

^  ^  °°^  j  pOO 

P{r,z)  =  Y.Y,  K’J (r’ 2)  -  P‘J (r> *)  =  9  /  Gij (z;p)E(0l) (wpr  )w'2p dp  .  (3) 

:=1  j= 0  i  J-'x> 

It  is  tempting  to  regard  the  PUj  as  ’generalized  rays’  and  consider  the  decomposition  (3)  as  an  exact 
extension  of  ray  theory  to  finite  frequencies.  The  traveling-wave  expansion  has  typically  been  utilized  in 


this  way  (e.g.,  [5]),  and  this  works  well  in  the  propagating  regime  for  wliich  the  WKBJ  approximation 
[6,  Box  9.3]  in  conjunction  with  the  stationary-phase  and  related  approximations  [7,  Appendix]  indicate 
dominant  contributions  from  the  neighborhoods  of  certain  well-defined  p  values.  Cf.  classical  ray  theory 
with  Airy-function  corrections  at  smooth  caustics  [8,  sec.  9-4], 

Here,  we  consider  the  traveling-wave  expansion  for  modeling  low-level  creeping-wave  diffraction  fields, 
which  turns  out  to  be  much  more  delicate,  hi  [9],  we  elucidated  the  importance  of  the  evanescent 
regime  of  wavenumbers  in  this  context.  It  turns  out,  however,  that  the  integrals  (over  p)  that  appear 
in  (3)  do  not  in  general  converge  at  infinity.  Hence,  the  decomposition  is  not  even  well  defined! 

As  it  stands,  the  decomposition  (3)  is  thus  not  useful  for  modeling  grazing-ray  diffraction  beyond 
boundary-produced  caustics.  For  single  grazing-ray  diffraction,  with  only  one  interaction  with  the 
boundary,  Pierce  [8,  sec.  9-5]  has  described  how  the  shadow-zone  field  can  be  represented  as  a  residues 
series  w'hose  leading  term  can  be  interpreted  physically  as  ray  shedding  by  a  creeping  wave.  Application 
of  this  residues  series  to  aeroacoustics  can  be  found  in  [10],  for  example. 

In  [11]  and  [12],  and  references  therein,  we  proposed  recombinant  (or  reduced)  traveling-wave  expan¬ 
sions  which  can  be  used  for  exact  modeling  of  single  as  well  as  multiple  grazing-ray  diffraction.  The 
convergence  criteria  are  relaxed,  allowing  specification-depth  location  at  the  point  of  minimum  sound 
speed  in  the  medium,  by  which  trapped  modes  are  avoided.  Physically,  ray  contributions  arise  from 
points  p  of  stationary  phase  while  ray  shedding  in  the  shadow'  zone  arises  from  leaky  modes. 

We  limit  ourselves  to  an  example  with  a  sound  channel  according  to  Fig.  la  with  source  and  receiver 
depths  at  10  m  and  50  m,  respectively.  The  subbottom  below  80  m  is  a  homogeneous  fluid  half-space 
with  density  2  kg/dm3,  velocity  2.0  km/s  and  absorption  0.7  dB/wavelength.  Transmission  loss  curves 
(re  1  m)  for  component  fields  PV2,j  are  shown  in  Fig.  lb.  The  recombinant-type  quantity  Pn,j  agrees 
with  the  sum  of  Pij  and  P-i.j  whenever  the  latter  are  well  defined. 


Figure  1:  a)  Sound-speed  profile,  b)  Recombinant-type  field  components  Pi2j  at  3  kHz.  • 

Just  as  for  the  downward  refracting  example  studied  in  [12],  the  recombinant-type  field  decomposition 
does  a  good  job  of  isolating  physical  ray  energy  to  the  term  with  the  appropriate  cycle  number  j.  The 
grazing  rays  that  start  upwards  from  the  source  and  reach  the  receiver  from  above  (type  i=2  according 
to  [5])  get  the  range  of  0.613,  2.119,  and  3.624  km  for  the  cycle  numbers  j=0,  1,  and  2,  respectively. 
These  ranges  are  clearly  related  to  the  shadow-zone  ranges  in  Fig.  lb,  beyond  wliich  the  pertinent 
P 12, j  components  decrease  steadily. 

Individual  leaky-mode  contributions  for  our  sound-channel  medium  example  are  studied  in  Fig.  2,  and 
also  in  Fig.  3a.  In  this  case,  it  is  mode  21  (as  numbered  from  the  right  in  Fig.  2a)  that  is  closest 
to  the  slowness  1/1.5  s/km  of  the  surface-grazing  ray,  but  the  imaginary  mode  slowness  parts  have 
an  increasing  trend  already  from  mode  1.  Although  mode  1  must  dominate  the  very  weak  P12,o 
field  at  very  long  ranges,  modes  close  to  the  surface-grazing  ray  slowness  w'ill  indeed  be  important  at 
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moderate  ranges.  In  Fig.  2b,  mode  20  provides  the  largest  contribution  at  r  =  0.3  krn,  whereas  mode 
17  dominates  at  r  =  0.6,  0.9,  and  1.2  km.  The  results  in  Fig.  3a  confirm  the  suitability,  for  the  purpose 
of  computational  efficiency,  of  classical  ray  theory  at  close  ranges  and  representation  with  a  limited 
number  of  leaky  modes  at  a  shadow  boundary  and  beyond. 


a)  s/kiri 


0.005- 


I _ , . . . •  0.0001 

0.66  0.67  s/km 

Figure  2:  a)  Mode  slownesses  in  the  complex  slowness  plane  relevant  to  a  leaky-mode  represen¬ 
tation  of  the  Fi2,o  field  in  Fig.  lb.  b)  Contributions  of  individual  modes  to  the  corresponding 
Pi2,o  field  at  four  different  ranges  r:  0.3,  0.6,  0.9,  and  1.2  km.  The  horizontal  axes  are  for 
mode  number,  according  to  decreasing  real  part  of  horizontal  slowness.  • 

For  components  with  cycle  number  j  >  1,  and  also  some  with  j  =  1,  the  additional  singularities  from 
7/1  and  cause  higher-order  poles.  The  corresponding  residues  will  involve  higher-order  derivatives 
that  do  not  seem  to  be  easily  computed.  To  be  able  to  conveniently  handle  all  field  components  under 
these  circumstances,  those  with  cycle  number  j  >  1  in  particular,  we  may  use  wavenumber  integration 
along  short  paths  in  the  complex  p  plane  to  capture  the  crucial  leaky  inodes. 

For  our  sound-channel  example,  the  corresponding  P 1  ‘2 results  arc  shown  in  Fig.  3b-d,  along  with  exact 
results  and  classical  ray-theory  results.  Again,  it  is  apparent  that  ray  theory  is  useful  at  close  ranges, 
whereas  leaky-mode  Green’s  function  expansion  is  appropriate  at  shadow  boundaries  and  beyond. 
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Figure  3:  a)  Transmission  loss  (re  1  m)  for  Pi2)0  as  in  Fig.  lb,  but  only  including  certain  leaky 
modes  from  Fig.  2.  Three  mode  curves  are  given,  with  increasingly  improved  agreement  to 
the  complete  P12]o  field.  They  were  obtained  as  coherent  sums  of  the  following  mode  fields: 
mode  17,  modes  12-21,  modes  1-100  in  Fig.  2b.  Two  additional  curves  are  given:  the  ’exact’ 
result  from  Fig.  lb  (displaced  40  dB  downwards  for  clarity)  and  the  classical  ray-theory  result 
(displaced  80  dB  downwards).  b)-d)  Transmission  loss  (re  1  m)  for  Pi2;0  (upper  right),  Pi2,i 
(lower  left),  Pi2,2  (lower  right),  as  obtained  by  wavenumber  integration  along  a  short  path  to 
enclose  the  horizontal  slownesses  of  the  crucial  leaky  modes.  The  corresponding  ’exact’  result 
from  Fig.  lb  (displaced  40  dB  downwards)  and  the  classical  ray-theory  result  (displaced  80 
dB  downwards)  are  also  shown  in  each  case.  • 
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INTRODUCTION. 


The  aim  of  the  present  work  is  to  demonstrate  the  practical  applications  of  a  high-power  ultrasound 
transducers  mounted  on  a  filtering  spiral  module.  This  construction  has  been  tested  in  a  vine 
production  plant  for  preventing  the  blocking  of  microfilters  by  organic  materials  and  for  increasing 
the  flow  rate  through  the  filters  and  the  time  between  the  filter  cleanings.  A  numerical,  FE  and  BE 
based  model  for  calculation  of  the  resonance  vibrations  of  a  liquid  filled  spiral  module  formed  the 
basis  for  the  design  of  the  transducer  attachment  construction.  During  the  laboratory  experiments 
the  resonance  response  of  the  filter  module  was  measured  carefully  and  compared  with  theoretical 
predictions.  The  properties  of  the  filtering  spiral  module  and  the  results  of  the  tests  are  discussed  in 
details.  This  work  was  financed  through  the  EU  Project  WAMBIO  PL96-3257  (FAIR  Programme). 

BACKGROUND  FOR  THE  RESEARCH. 


In  vine  filtering  processes  a  spiral  filtering  module  is  commonly  used.  In  order  to  find  the  optimum 
solution  for  design  of  the  attachment  of  ultrasound  transducers  to  the  filtering  module  careful  choice 
of  the  attachment  point  should  be  made.  The  aim  of  such  a  construction  is  to  excite  the  natural 
vibration  modes  of  the  filter  housing  (which  is  shaped  as  a  cylinder  with  closed  ends)  as  a  whole, 
where  flexural  vibrations  can  travel  along  the  cylindrical  tube. 

THEORETICAL  MODEL 

In  our  case  the  simplest  model  for  the  spiral  wound  filter  is  the  cylindrical  tube  with  two  closed 
ends.  In  agreement  with  the  technical  specification  the  length  of  the  steel  tube  is  960mm,  with  a 
inner  diameter  of  104  mm  and  a  wall  thickness  of  h=2,5  mm.  Tube  mean  radius  is  53,25  mm 
Young  modulus  E  of  the  tube  material  is  20-1010  N/m2;  the  shear  modulus  G  is  7.69-1010  N/m2  and 
the  Poisson’s  ratio  v  =  0.245.  Figure  1  shows  a  sketch  of  the  flexural  waves  travelling  along  the 
cylindrical  tube. 
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Figure  1.  Flexural  wave  propagation  on  a  thin-walled  cylindrical  tube. 


Two  important  parameters  for  the  definition  of  the  resonance  properties  of  the  cylindrical  tube  are: 

•  The  nondimentional  frequency  Q=ooa/ci,  where  to  is  the  angular  frequency,  a  is  the  mean 
radius  of  the  tube  and  ci  is  the  longitudinal  wave  speed. 

•  The  nondimentional  thickness  parameter  (3=h/(Vl2)a 

The  first  parameter  indicates  the  frequency  at  which  an  antisymmetric  resonance  mode  can  occur 
(n  is  the  mode  number) 

If  we  suppose,  that  the  thickness  of  the  tube  wall  is  small  enough  compared  to  the  mean  radius,  (i.e. 
h  «  a),  then  the  frequency  equation  for  the  antisymmetric  resonance  modes  of  the  cylindrical  tube 
will  be: 

Q2n=|32n4[l-0.5  (l/(l-v){(4-v)/n2-  (2+v)/n4})]. 

Simple  calculations  will  give  us  all  resonance  frequencies  corresponding  to  the  different  mode 


numbers. 

For 

n=5 

Qn=  (3-23,76=0.0751 

oo=8808  and 

f=1402 Hz 

For 

n=10 

Qn=  (3-97,53=0.3082 

co=36155  and 

f=5754  Hz 

For 

n=19 

f =21 223  Hz 

For 

n=20 

f=  23525  Hz 

For 

n=21 

f=25945  Hz 

For 

n=22 

f=28483  Hz 

For 

n=23 

f =31 140  Hz . 

Note,  that  our  model  tube  is  empty  (if  filled  with  water  the  resonance  frequencies  will  be  reduced). 

If  a  flange  is  introduced  (as  it  is  in  the  case  of  a  spiral  filter),  the  resonance  frequencies  will  be 
increased. 

EXPERIMENTAL  TESTS 

From  the  theoretical  investigations  of  the  vibration  modes  of  the  simple  model  of  a  steel  tube,  it  was 
concluded  that  the  excitation  of  the  steel  tube  as  a  whole  (also  comprising  the  spiral  filter 
module)  can  potentially  be  used  for  obtaining  the  necessary  high  levels  of  ultrasound  intensity 
inside  the  filter  unit  and  that  a  frequency  around  25-27  kHz  could  be  used. 
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A  94  cm  x  10  cm  filter  module  containing  a  standard  plastic  microfilter  membrane  MPPS-U002 
was  assembled,  filled  with  water  and  tested  in  the  laboratory.  A  RESON  standard  TC  5003 
transducer  was  used  to  excite  the  tube  vibrations.  An  approx.  25  kHz  continuous  sinusoidal  signal 
was  generated  by  a  TCE  7702  signal  generator  and  the  signal  was  amplified  in  a  HVPA-01  power 
amplifier. 

The  flexural  vibrations  were  picked  up  by  a  B&K  accelerometer  and  monitored  by  a  HP  54602A 
oscilloscope. 

The  accelerometer  was  fixed  to  the  outside  tube  wall  using  bee  wax,  and  a  working  line  for 
transducer  excitation  was  drawn  from  the  accelerometer  to  the  end  of  the  tube. 

The  head  part  of  the  TC  5003  transducer  had  a  diameter  of  60  mm,  while  the  local  resonance  points 
of  the  tube  was  expected  to  have  a  size  of  only  5  mm.  Therefore,  the  ultrasonic  vibrations  were 
applied  to  the  steel  tube  by  turning  the  transducer  around  and  by  pressing  the  transducers  bolt  head 
vertically  against  the  tube  wall. 

Registration  of  the  tube  wall  vibration  modes 

Each  local  maximum  response  location  was  recorded.  The  schematic  picture  below  shows  the 
experimental  set-up  used  during  the  measurement  of  the  resonance  response  along  the  tube  module 
with  the  indication  of  the  position  of  the  transducer  and  the  positions  of  the  recorded  maximum  of 
the  response  signals.  The  dimensions  correspond  to  a  working  frequency  of  27  KHz 
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During  the  tests  the  accelerometer  was  fixed  in  five  different  positions  on  the  inlet  /outlet  side  of  the 
tube,  and  the  whole  measurement  procedure  was  repeated  for  each  accelerometer  position. 

The  procedure  was  repeated  for  3  different  frequencies:  24  kHz,  25  kHz  and  27  kHz  with  both 
transducer  and  accelerometer  operating  at  the  opposite  side  of  the  inlet/outlet  side  of  the  tube. 


DISCUSSIONS  OF  THE  EXPERIMENTAL  RESULTS  AND  THE  CONCLUSIONS 

Careful  analysis  of  experimental  data  allowed  us  to  draw  several  conclusions: 

1.  The  best  performance  was  obtained  using  a  27  kHz  signal  to  excite  the  filtering  module.  The 
tube  response  signal  detected  by  accelerometer  have  had  a  more  or  less  periodic  structure  with 
the  repetition  rate  of  the  maximum  value  of  the  response  signal  after  every  5  cm.  The 
amplitude  level  at  maximum  response  was  about  the  same  along  the  whole  tube,  around  240- 
280  mV.  It  means  that  it  was  possible  to  excite  the  vibration  mode  of  the  full  spiral  filtering 
module  by  propagation  of  flexural  wave  signals  along  the  module. 

2.  The  tested  frequency  of  25  kHz  shows  a  periodical  structure  for  the  maximum  signal  response 
in  the  axial  direction  of  the  tube  with  a  repetition  rate  of  approx.  4  cm  between  the  maximum 
values.  Maximum  amplitude  for  the  response  signal  was  recorded  to  140  mV  (only  half  of  the 
amplitude  measured  at  27  kHz). 

3.  For  the  24  kHz  signal  no  periodic  behaviour  was  founded.  A  maximum  response  accelerometer 
signal  of  100  mV  was  recorded. 

4.  Calculation  of  the  resonance  frequencies  using  the  simple  model  for  the  spiral  wound  filter 
shown  earlier  gives  the  value  for  the  eigen-frequency  of  the  empty  tube  to  be  28483  Hz.  If  we 
take  into  account  that  the  simple  model  can  only  give  an  indication  of  the  magnitude  of  the 
actual  flexural  wave  frequency,  the  experimental  value  of  the  best  performance  frequency  of 
27  kHz,  is  very  close  to  the  theoretically  predicted  one.  The  reduction  from  the  calculated 
empty  tube  frequency  of  28.48  kHz  to  the  measured  27  kHz  for  optimal  operation  should  be 
expected  due  to  the  presence  of  liquid  and  filter  module  in  the  tube  during  the  experiments. 
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5.  Efficient  transducer  attachments  to  the  spiral  wound  filter  housing  for  the  excitation  of  the 
flexural  vibration  modes  of  the  tube  wall  thus  leading  to  an  ultrasonic  cleaning  effect  on  the 
filter  module  can  be  achieved. 

6.  Several  technical  problems  should,  however,  be  solved  before  a  successful  application  of  the 
ultrasound  transducers  to  the  cleaning  of  the  filter  module  can  be  carried  out: 

•  New  construction  for  the  attachment  of  the  transducer  to  the  cylindrical  filter  surface  should 
be  considered  for  maximum  transfer  of  power  to  the  tube  wall. 

•  Optimisation  of  the  transducer  performance  should  be  considered  as  an  alternative  as 
calculations  have  shown,  that  the  use  of  lower  frequencies  could  lead  to  a  more  efficient 
excitation  of  the  tube  wall  vibration. 

•  Losses  at  high  frequencies,  which  are  about  90%,  may  be  reduced  by  using  low  frequency 
transducers  (fundamental  frequency  of  about  1,5  kHz)  and  by  improving  the  attachment 
procedure. 
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Sound  radiation  from  unbaffled  loudspeakers 
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Abstract 

The  near-field  of  an  unbaffled  loudspeaker  is  studied  A  model  employing  edge  diffraction  is  used,  and  this  gives 
an  analytic  expression  for  the  impulse  response  of  a  circular  piston  taking  the  direct  sound  and  the  first-order 
diffraction  into  account.  Higher-order  diffraction  is  approximated  by  a  rest  function.  Results  show  that,  as 
expected,  an  unbaffled  loudspeaker  element  gives  a  less  homogeneous  sound  field  for  movements  along  the  sym¬ 
metry  axis  than  a  baffled  loudspeaker  element  does.  On  the  other  hand,  the  unbaffled  loudspeaker  also  leads  to  a 
lower  crosstalk  between  the  two  ears  of  a  listening  person  sitting  very  close  to  the  loudspeaker  element  than  a 
baffled  loudspeaker  does. 

Introduction 

Loudspeakers  are  usually  placed  in  baffles  or  enclosures  in  order  to  transform  their  inherent  dipolar  operation 
into  a  monopolar  operation  with  its  higher  efficiency.  For  some  applications,  however,  unbaffled  loudspeakers 
might  be  used.  A  recently  presented  active  noise  control  system  for  the  creation  of  a  so-called  Silent  Zone 
around  a  person's  head,  see  Fig.  1  [1],  uses  unbaffled  loudspeakers.  For  such  applications  the  nearfield  of  the 
unbaffled  loudspeaker  must  be  analyzed.  The  purpose  of  this  paper  is  to  study  using  a  model  of  a  vibrating,  free 
circular  piston.  An  impulse  response  method  using  the  edge  diffraction  concept  is  used  for  the  study  [2].  First, 
the  edge  diffraction  impulse  response  method  is  reviewed  briefly.  Then  the  case  of  a  piston  is  modelled  as  a 
superposition  of  point  sources  and  numerical  results  are  presented  for  this  case. 


Analytic  edge  sources 

The  scattering  from  an  object  with  rigid,  plane  surfaces  can  be  formulated  as  a  sum  of  geometrical  acoustics 
components,  i.e.,  the  direct  sound  and  specular  reflections,  and  edge  diffraction  components.  Edge  diffraction  is 
caused  by  the  waves  that  are  scattered  off  the  straight  or  curved  edges  of  the  finite  surfaces  and  can  be  represen¬ 
ted  by  secondary  sources  placed  along  the  edges  [2,3].  These  sources  radiate  spherically  in  addition  to  a 
directivity  function  which  depends  on  the  angle  of  the  incident  wave  from  the  point  source  towards  the  edge 
point  (angles  a  and  0sin  Fig.  2),  and  also  on  the  radiation  angle  (angles  yand  0Kin  Fig.  2).  The  decomposition 
can  be  expressed  as  a  line  integral  along  the  edge  of  the  wedge  so  that  the  sound  pressure  IR  is 
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where  c  is  the  speed  of  sound,  v  is  the  so-called  wedge  index  which  equals  jt/0w,  and  m  and  /  are  the  distances  to 
and  from  a  point  z  along  the  edge  as  indicated  in  Fig.  2.  The  factor  (3  is  a  directivity  function  which  has  the  form 


sin 

V\JZ  ±  % 

+eR) 

4 

±+cosh 

f  ,  _i  1-sinasiniO 

vcosh  - - - 

cosacosy  J 

-cos 

v(jr±e5+0fl)] 

(2) 


The  summation  includes  four  terms,  with  all  possible  combinations  of  the  two  signs  indicated  in  the  sine  and 
cosine  factors,  and  with  the  sign  combinations  in  the  numerator  and  denominator  being  the  same.  It  is  clear  that 
P  can  be  interpreted  as  a  directivity  factor  since  it  includes  only  the  incidence  angles  a  and  0S,  and  the  reradia¬ 
tion  angles  y  and  0R.  For  an  IR  as  in  Eq.  (1),  the  source  signal  quantity  is  p0v4(r)/(4ji)  where  p0  is  the  density  of 
the  fluid  and  A(t)  is  the  volume  acceleration  of  the  point  source.  This  choice  of  source  signal  gives  a  free-field  IR 
which  is  /i  (t)  =  b(t-R/c)/R. 

The  expressions  in  Eqs.  (1)  and  (2)  can  be  extended  to  curved  edges  as  long  as  the  integration  variable  z  is 
interpreted  as  running  along  the  curved  edge.  Also,  finite  wedges  can  be  handled  directly  by  carrying  out  the  in¬ 
tegration  over  the  range  of  z  which  corresponds  to  the  finite  wedge.  Furthermore,  the  extension  to  higher  orders 
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FIG.  1  Illustration  of  an  active  noise  control  system 
implementing  a  silent  zone  around  a  person's  head  [1], 
The  loudspeaker  elements  are  unbaffled. 
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FIG.  2  Illustration  of  the  various  parameters  for  the 
diffraction  from  a  point  source  S  via  a  point  along  the 
edge  of  a  wedge  to  a  receiver  point  R.  The  two  fictive 
planes  created  by  the  edge  and  the  source  and  the 
receiver,  respectively,  are  indicated  by  the  dashed  lines. 


of  edge  diffraction,  which  might  be  needed  to  include  in  many  situations,  is  relatively  straightforward,  see  [2],  It 
should  be  noted  that  wedges  which  have  one  or  two  pressure-release  planes,  rather  than  two  rigid  planes,  can  be 
handled  simply  by  changing  the  sign  of  one  or  more  of  the  four  terms  in  the  sum  in  Eq.  (2).  Having  access  to 
analytical  directivity  functions  for  the  secondary  edge  sources  offers  the  possibility  to  derive  analytical  expres¬ 
sions  for  the  edge  diffraction  IR  if  the  angles  in  Eq.  (2)  can  be  expressed  in  terms  of  the  integration  variable  z. 
This  is  demonstrated  in  the  next  Section  for  the  case  of  a  circular  piston.  Finally,  the  expression  in  Eq.  (1)  is  very 
well  suited  for  numerical  implementation  since  the  integrand  is  well  behaved.  Consequently,  first-order 
diffraction  can  be  calculated  numerically  for  arbitrary  geometries  with  high  accuracy  and  very  little  computa¬ 
tional  effort.  It  should  be  noted  that  this  is  exactly  what  the  so-called  Wedge  Assemblage  (WA)  method  does* 
[4],  The  WA  method  is  based  on  Medwin's  interpretation  of  the  Biot-Tolstoy  method,  [5],  and  Medwin's  method 
in  [3]  is,  numerically,  identical  with  the  presented  method  for  first-order  diffraction  but  is  approximate  for 
higher-order  diffraction. 


The  unbaffled  piston 

A  vibrating  piston  can  be  modelled  as  a  rigid  thin  disc  and  a  continuous  distribution  of  point  sources  on  both 
sides  of  the  scattering  disc,  with  opposite  polarity.  The  total  sound  field  can  then  be  described  by  an  impulse 
response  which  is  a  sum  of  the  direct  sound  IR  (from  the  frontal  side  of  the  disc),  and  the  edge  diffraction  IR:s 
for  increasing  orders.  This  is  a  special  case  with  v  =  1/2  and  6S  =  0.  For  that  case,  illustrated  in  Fig.  3,  the 
directivity  function  can  be  simplified  into 

p  ScoshMcos-^  8V2cosh(v,?),/IT^ 

cosh 2  (vtj)  -  sin 2  2cosh2(v7/)-l  +  cos0/?  ^ 


where 


cosh(vrj)  = 
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In  Eq.  (3),  the  effect  of  the  sources  on  the  rear  side  of  the  disc,  with  opposite  polarity,  have  been  taken  into 
account  since  their  contribution  is  exactly  the  same  as  the  sources  on  the  frontal  side.  For  a  continuous  distribu¬ 
tion  of  sources,  the  edge  diffraction  component  can  be  derived  by  integrating  the  contributions  from  all  source 
positions, 

htot,D\(t)  =  (5) 

s 


where  S  is  the  piston's  area,  and  /iDI(f)  is  the  first-order  diffraction  IR  for  each  source  position.  If  we  study  only 
on-axis  receiver  positions,  this  integral  can  be  solved  analytically  and  the  result  is 


where  the  receiver  is  at  a  height  z,  which  gives  /  =  Va2  +  z2  ,  a  being  the  radius  of  the  piston.  The  function  w(t) 
is  a  window  function  which  is  zero  before  the  first  possible  arrival  of  sound  via  an  edge  point,  and  zero  after  the 
last  possible  arrival  of  sound.  Between  those  two  points  in  time,  w(r)  =1. 


108 


R 


/• 
/  • 


FIG.  3  A  circular  disc  with  a  source  point  S  on  the 
disc's  surface  and  a  receiver  point  off-axis. 
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FIG.  5  The  transfer  function  of  the  IR  in  Fig.  4,  but 
with  half  the  positive  pulse  in  Fig.  4  removed. 
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FIG.  4  IR  for  an  unbaffled  piston  of  radius  0.1  m  and 
a  receiver  on-axis  (i.e.,  with  y  =  0)  at  a  height  of  0.1  m. 
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FIG.  6  The  rest  function,  which  replaces  second-  and 
higher-order  diffraction. 


Fig.  4  illustrates  the  IR  which  includes  the  direct  sound  of  the  piston,  and  the  edge  diffraction  IR  according  to 
eq.  (6).  Second-  and  higher  diffraction  components  can  not  be  derived  analytically  so  numerical  integration  is 
necessary.  Here,  a  simplified  approach  has  been  chosen  where  the  second-  and  higher-order  diffraction  is  re¬ 
placed  by  an  approximate  "rest  function".  A  reasonable  estimate  of  this  rest  function  can  be  found  since  we 
know  that  the  total  impulse  response  from  the  piston  must  asymptotically  go  towards  free-field  radiation  for  very 
low  frequencies.  In  Fig.  4,  the  first  positive  rectangular  pulse  is  the  baffled  piston  radiation,  given  by  radiation 
into  2n,  and  thus  with  twice  the  amplitude  of  free-field  radiation.  The  sum  of  half  the  first  rectangular  pulse,  plus 
all  orders  of  diffraction  must  then  be  a  high-pass  filter.  Fig.  5  illustrates  this.  Stated  another  way,  the  scattering 
of  a  rigid  disc  must  be  a  high-pass  filter,  and  here  we  have  a  scattering  case  with  a  free-field  radiating  distribu¬ 
tion  of  point  source  immediate  close  to  the  rigid  disc.  We  can  choose  a  rest  function  which  is  shaped  as  in  Fig.  6 
since  this  shape  should  match  the  second-order  diffraction  quite  closely.  It  is  clear  at  which  times  the  second- 
order  diffraction  must  start  and  stop-  The  area  of  the  rest  function  is  chosen  so  that  the  total  response  of  reflec¬ 
tion  +  first-order  diffraction  (eq.(6))  and  the  rest  function  is  0  at  0  Hz.  This  is  clearly  a  very  simplified  approach 
but  it  can  give  very  useful  estimates  of  the  transfer  function  also  at  low  frequencies  and  uses  explicit  simple 
expressions. 


Numerical  examples 

For  an  application  such  as  in  Fig.l,  three  aspects  of  the  near-field  are  especially  interesting:  how  homogeneous 
the  sound  field  is  if  the  listener  moves  the  head,  how  strong  the  cross-talk  is  (from  the  nearest  to  the  furthest  ear) 
and  how  efficient  an  unbaffled  piston  is  compared  to  a  baffled  piston.  These  points  were  studied  in  the  free-field 
in  front  of  the  unbaffled  piston,  i.e.,  the  reflecting/scattering  of  the  head  was  ignored.  At  low  frequencies  this 
assumption  is  quite  acceptable. 

The  homogeneity  in  front  of  the  piston  is  important  in  an  active  noise  control  system  since  such  a  system 
aims  at  supplying  an  out-of-phase  sound  wave  with  a  specific  amplitude  and  phase  controlled  as  accurately  as 
possible.  This  can  be  accomplished  only  in  a  single  point  but  if  the  sound  field  is  homogeneous,  the  con¬ 
trolled/silent  zone  will  be  larger.  Fig.  7  shows  the  transfer  function  (TF)  magnitude  for  the  frequency  of  94  Hz, 
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FIG.  7  The  TF  magnitude  at  94  Hz  for  a  circular  FIG.  8  The  TF  magnitude  difference  for  a  baffled 

piston  of  20  cm  diameter  with  the  receiver  on-axis.  piston  and  with  two  on-axis  receiver  positions  20  cm 

apart.  Parameter  is  the  distance  between  the  piston  and 
the  closest  microphone. 


Frequency  [Hz] 

FIG.  9  Same  as  in  Fig.  8  but  for  an  unbaffled  piston. 

as  function  of  distance  from  the  piston,  for  receiver  positions  on-axis.  Results  are  shown  for  both  a  baffled  and 
an  unbaffled  piston.  As  can  be  seen,  a  baffled  piston  does  actually  give  a  sound  field  which  is  varying  less  for 
distances  around  10-20  cm  from  the  piston.  It  is  also  clear  that  the  baffled  piston  gives  around  20  dB  more  output 
than  the  unbaffled  piston.  This  is  an  important  issue  as  well  since  the  loudspeaker  element  might  need  to  produce 
high  levels  of  noise  cancelling  sound. 

An  active  noise  control  system  might  use  two  independent  channels  or  use  some  kind  of  coupled  control  of 
the  two  silent  zones  around  the  two  ears.  The  crosstalk  is  an  all-important  parameter  for  this;  if  one  loudspeaker 
gives  a  strong  contribution  both  to  its  own  ear  and  to  the  ear  across  the  head,  the  system  might  have  to  use  a 
coupled  control.  On  the  other  hand,  if  the  crosstalk  is  very  low,  independent  channels  should  be  perfectly  ade¬ 
quate.  Figs.  8  and  9  show  the  magnitude  of  the  crosstalk  for  the  baffled  and  the  unbaffled  piston,  respectively. 
The  crosstalk  was  calculated  simply  by  comparing  the  TF  magnitude  for  two  receiver  positions  20  cm  apart,  on- 
axis.  This  corresponds  approximately  to  the  distance  between  the  ears  if  one  takes  into  account  the  longer  sound 
path  around  the  head.  In  reality,  the  crosstalk  might  be  different  than  this  since  the  head  does  have  some,  if 
weak,  shadowing  effect.  At  the  same  time,  the  other  loudspeaker  might  act  as  a  reflector  and  thus  increase  the 
crosstalk  magnitude.  Ignoring  these  effects  as  a  first  approximation,  it  can  be  seen  that  the  crosstalk  for  the  baff¬ 
led  piston,  in  Fig.  8,  is  quite  much  stronger  than  for  the  unbaffled  piston,  in  Fig.  9.  This  is  especially  true  if  the 
closest  ear  is  very  close  to  the  loudspeaker.  At  high  frequencies,  the  unbaffled  and  the  baffled  piston  behave 
identically. 

The  crosstalk  also  has  important  consequences  in  binaural  sound  reproduction  over  loudspeakers.  So-called 
crosstalk  cancellation  filters  are  needed  if  the  crosstalk  is  strong,  and  such  filters  are  quite  sensitive  to  the 
listening  person's  head  size  and  positioning  of  the  head.  Consequently,  a  reproduction  system  using  closely 
mounted  unbaffled  loudspeakers  would  not  need  any  cross-talk  cancellation  filter. 

Conclusions 

The  near-field  of  an  unbaffled  loudspeaker  has  been  studied  and  compared  with  that  of  a  baffled  loudspeaker.  An 
edge  diffraction  model  was  used  for  the  numerical  examples.  It  was  shown  that,  as  expected,  the  dipolar  nature 
of  the  unbaffled  loudspeaker  element  gives  a  less  homogeneous  nearfield,  but  with  a  sigificantly  lower  crosstalk 
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between  two  ear  positions  than  a  baffled  loudspeaker  would  have  given.  The  output  of  the  unbaffled  loudspeaker 
might  be  20  dB  lower  than  the  baffled  loudspeaker  around  100  Hz,  or  14  dB  lower  than  a  loudspeaker  in  a  small 
enclosure. 

The  calculation  model  used  gave  an  analytic  expression  for  the  first-order  diffraction  and  an  approximate  rest 
function  replacing  the  higher  orders  lead  to  a  very  efficient  numerical  method. 
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